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ABSTRACT 


Research  to-date  towards  an  understanding  of  human  biped 
locomotion  has  been  primarily  experimental  in  nature,  largely  due  to 
the  complexity  of  the  process.  In  view  of  the  new,  exciting  possibilities 
of  programmed  electro-stimulation  of  paralyzed  extremities  to  restore 
locomotion,  a  critical  study  at  the  theoretical  level  is  greatly  warranted 
Optimal  programming  and  modern  control  theory  offer  a  new  approach 
to  the  study.  First,  it  is  proposed  that  normal  walking  obeys  a 
certain  "principle  of  optimality Next,  at  the  dynamic  level,  modern 
control  theory  is  used  to  derive  the  optimal  moment  profiles  which 
actuate  the  locomotor  elements  to  synthesize  the  observed  patterns  of 
the  normal  gait. 

Development  of  the  problem  structure  relies  closely  on  the  func¬ 
tional  characteristics  of  the  biped  gait,  particularly  the  ideas  of  distinct 
phasic  activities  and  the  associated  temporal  patterns  of  a  walking  cycle 
The  result  is  a  multi-arc  programming  problem  with  three  stages. 

Each  stage  involves  dynamic  constraints  which  reflect  the  particular 
nature  of  the  phasic  activity.  Activity  in  the  stance  phase  is  described 
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by  equality  constraints  on  the  "states”  while  the  swing  phase  is  charac¬ 
terized  by  inequality  state  constraints.  A  novelty  of  the  approach  is 
that  the  theory  could  be  used  to  study  walking  behavior  under  different 
environmental  conditions,  such  as  walking  up-stairs  or  over  a  hole. 
Joining  of  the  arcs  is  arranged  in  such  a  way  as  to  maintain  continuity 
of  certain  trajectories  as  well  as  repeatability  of  motion. 

A  distinct  feature  of  the  present  approach  which  differs  from 
other  studies  is  the  presence  of  a  minimizing  performance  criterion. 
Based  on  external  characteristics  of  muscles  and  certain  assumptions 
regarding  normal  locomotion,  a  simple  quadratic  type  of  performance 
index  is  proposed.  This  performance  criterion  is  meaningful  in  that  it 
is  shown  to  be  proportional  to  the  mechanical  work  done  during  normal 
walking. 

Invoking  the  necessary  conditions  of  optimal  control  theory,  a 
multi -point  boundary  value  problem  is  obtained.  A  penalty  function 
technique  is  then  employed  for  iterative  numerical  solution  on  a 
digital  computer.  Using  the  parameters  available  in  the  literature, 
simulation  results  are  obtained  which  agree  well  with  the  experimental 
studies  performed  by  Eberhart  and  his  associates.  Furthermore, 
certain  refined  features  are  obtained  which  are  not  available  in 
previous  studies.  Success  in  applying  optimal  programming  techniques 
to  human  locomotion  could  yield  better  design  procedures  for  prostheses 
and  could  allow  eventual  realization  of  the  dream  of  programmed 
stimulation  of  many  paralyzed  persons. 
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A.  Preliminary  Considerations 


1 .  Motivation. 

Human  locomotion  displays  a  high  degree  of  complexity,  organi¬ 
zation,  and  efficiency.  The  versatility  and  "gracefulness”  of  the  motion 
has  fascinated  man  for  a  long  time.  According  to  Bernstein  [1], 
Leonardo  da  Vinci  was  probably  the  first  to  make  scientific  observations 
about  human  locomotion  although  real  quantitative  studies  are  only  of 
recent  vintage.  In  addition  to  scientific  curiosity,  the  phylogenetic 
development  of  the  biped  gait  from  the  push-pull  mechanism  of  the 
quadruped  is  also  of  interest  to  the  physical  anthropologist.  How  much 
better  is  a  biped  adapted  to  its  living  patterns  than,  say,  its  quadruped 
counterpart  To  answer  this  question  adequately,  one  has  to  explore 
the  functional  characteristics  of  the  biped  gait. 

An  important,  practical  goal  of  locomotion  study  is  a  better 
design  of  prostheses  for  disabled  persons.  In  this  and  other  areas 
dealing  with  human  limbs  and  their  substitutes,  a  common  consensus  is 
that  the  end  products  should  duplicate  the  performance  of  their  normal 
human  counterparts  as  closely  as  possible.  Thus,  quantitative 
knowledge  of  the  normal  biped  gait  serves  as  a  common  denominator 
for  performance  evaluation  and  improvement  of  the  various  innovations. 

The  modern  study  of  human  locomotion  began  in  the  1930’s  and 
has  drawn  greater  attention  since  then.  Significant  contributions  during 
this  time  include  the  works  of  Fenn  [2,  3],  Steindler  [4],  Elftman  [  5 - 8 ] , 
Liberson  [9],  Wilson  [10],  Close  [11],  and  the  Ebe rhart -Inman  group  at 
Berkeley  [12,  13].  Largely  due  to  the  complexity  of  the  subject,  results 
obtained  so  far  have  been  experimental  in  nature.  However,  in  view  of 
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the  recent  exciting  possibilities  of  restoring  the  locomotion  ability  of 
certain  paralyzed  limbs  through  programmed  electro  -  stimulation  [14- 
20],  a  critical  study  of  biped  locomotion  at  the  theoretical  level  is 
greatly  warranted.  After  all,  the  very  fundamental  question  of  M  what 
is  the  basic  mechanism  of  locomotion?”  or  ”what  are  the  laws 
governing  locomotion?”  cannot  be  answered  in  its  completeness  by 
experimental  study.  A  theory  which  attempts  to  embody  the  known 
facts  about  the  human  gait  would  be  an  important  step  in  this  field. 

Research  activities  to-date  appear  to  fall  into  two  broad  areas. 

The  first  area  is  concerned  with  the  functional  structure  and  operational 
behavior  of  the  locomotor  system  and  its  gait.  The  second  studies 
energy  expenditure,  oxygen  consumption  and  metabolism  in  general 
associated  with  walking.  These  two  areas  describe  the  important 
aspects  of  a  very  complicated  process.  With  the  tools  of  modern 
control  theory,  it  appears  possible  to  weld  these  two  areas  into  a 
’’functional”  theory  of  human  locomotion.  Such  a  study  leads  to  quanti¬ 
tative  investigation  and  useful  design  procedures  for  prosthetic  devices. 
To  this  end,  the  relevant  physiological  information  and  gait  features 
are  summarized.  This,  then,  sets  up  the  framework  of  our  theoretical 
study. 

2 .  Quantitative  Features  of  the  Biped  Gait 

a.  Physiological  Factors  [21] 

Locomotion  activity  and  other  voluntary  efforts  usually  involve 
decision  and  supervision  by  the  ’’higher  centers  ”  in  the  brain.  Fig.  1 
illustrates  a  typical  skeletal  muscle  and  its  neural  pathways.  The  skel¬ 
etal  muscle  is  the  actuating  unit  for  all  locomotion  activities.  The 
a  -motor  input  along  the  spinal  cord  represents  a  direct  input  from  the 
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FIG.  1  SKELETAL  MUSCLE  AND  NEURAL  PATHWAYS. 
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FIG.  2 


BLOCK  DIAGRAM  FOR  NEURAL- MUSCULAR  CONTROL  ACTION. 
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higher  centers  to  the  muscle  by  which  the  desired  action  can  be  initiated. 
It  is  believed  that  the  a  -motor  input  is  the  signal  used  for  voluntary 
actions.  A  voluntary  action  exists  in  the  conscious  thought  in  terms  of  a 
functional  rather  than  a  specific  description.  "Lower"  in  the  brain,  the 
concept  of  "voluntary  actions"  is  translated  into  a  programmed  sequence 
of  individual  muscle  actions,  based  on  the  person’s  experience  in  the 
past  in  performing  the  task.  The  "activation"  signal  travels  along  the 
a  -motor  pathway.  A  muscle  spindle  receptor  is  shown  inserted  in  a 
parallel  position  to  the  skeletal  muscle.  This  is  a  transduce r  device 
which  measures  the  length  of  the  muscle  and  its  sensitivity  is  adjusted 
by  the  upper  y-motor  neuron.  Results  of  the  measurement  are  fed 
back  to  the  thick  lower  a  -motor  neuron  at  the  synapse  where  further 
integration  of  activity  is  taken.  Thus  the  upper  a  -motor  neuron  as 
well  as  the  peripheral  sensory  neuron  of  the  muscle  receptor  control 
the  lower  a  -motor  neuron  and  subsequent  muscular  activity. 

The  simple  picture  of  neural  muscular  action  can  be  shown 
functionally  in  a  block  diagram  in  Fig.  2.  The  present  study  is  pri¬ 
marily  concerned  with  the  various  mechanical  moments  which  actuate 
the  locomotor  system  and  the  resulting  gait  patterns  it  produces.  By 
knowing  the  moment  histories,  one  can  then  concentrate  on  the  neural- 
muscular  transformation  (E.M.G.  activity  and  the  like)  for  the 
electrical  patterns  of  activity  and  eventually  the  neural  action.  This 
amounts  to  working  outwards  from  the  inner-most  loop  rather  than 
studying  the  entire  system  in  aggregate. 
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b.  The  Compass  Gait  [9,22] 

Under  the  action  of  the  actuating  moments,  the  motion  of  the 
locomotor  system  produces  the  patterns  of  the  biped  gait  which  can  be 
studied  experimentally  to  various  degrees  of  accuracy.  The  principal 
feature  underlying  the  biped  gait  is  the  compass  motion  of  the  two 
lower  extremities  about  alternate  centers  of  rotation.  Such  a  rotation 
of  the  body  about  alternate  centers  brings  about  an  overall  translatory 
motion.  In  Fig.  3,  suppose  that  the  left  leg  has  just  completed  its 
foreward  swing  and  come  into  the  restraining  portion  of  the  support 
(stance)  phase.  The  deploy  action  of  the  right  leg  pushes  the  body  for¬ 
ward  with  the  hip  describing  a  circular  arc  ABC  about  the  supporting 
left  foot.  At  passing  point  B,  the  center  of  gravity  of  the  body  is 
directly  above  the  femural  head  of  the  thigh.  After  that,  the  hip  pre- 
scribes  the  segment  BC  with  the  right  leg  in  swing  motion  under 
gravity.  At  C,  the  swinging  right  leg  restrains  the  falling  tendency 
and  goes  into  the  stance  phase  of  activity.  Concurrently,  the  left  leg 
in  the  late  stance  undergoes  deploy  and  continues  into  swing.  This 
results  in  a  second  arc  CDE  for  the  hip  with  the  role  of  the  left  and 
right  legs  interchanged.  When  this  is  completed,  the  entire  cycle  of 
the  "double  -  step  "  is  repeated.  The  compass  gait  results  in  a  series  of 
intersecting  circular  arcs  for  the  hip  trajectory  in  level  walking.  In 
actual  motion,  the  gait  is  made  more  graceful  and  efficient  during 
transfer  between  adjacent  arcs  by  additional  motion  of  the  legs  and 
pelvic  adjustment.  The  smoothed  action  results  in  an  almost  sinu¬ 
soidal  trajectory  for  the  hip. 


FIG.  3  THE  COMPASS  GAIT. 
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FIG.  4  TEMPORAL  PATTERNS  OF  BIPED  GAIT. 
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c.  Phasic  Activities  and  Temporal  Patterns  [4] 

The  smoothed  pattern  of  the  biped  gait  requires  a  good  deal  of 
coordination  in  the  activities  of  the  two  lower  extremities.  The  result 
is  that  there  are  two  characteristic  phases  of  activity  in  a  "double  -  step  n 
the  stance  (restraint,  support  and  deploy)  and  swing.  In  normal 
walking,  the  phasic  activities  of  one  leg  occur  with  regularity  and  are 
coordinated  with  those  of  the  other  leg.  Fig.  4  illustrates  such  a 
temporal  pattern.  Note  that  there  is  a  brief  double  stand  period 
(restraint/deploy)  when  both  legs  are  in  contact  with  the  ground.  In 
fast  walking,  approaching  running,  tne  double  stand  becomes  shorter. 

In  running,  there  is  even  a  brief  period  when  both  legs  are  in  the  air. 
The  mechanism  of  the  striding  mode  in  normal  walking  is  different 
from  that  of  the  running  mode.  Using  a  foot-switch  method,  it  has  been 
measured  that  the  stance  phase  occupies  roughly  60%  of  the  time 
required  for  a  double-step.  The  remaining  40%  is  swing  motion.  The 
deploy  (or  restraint)  portion  in  stance  takes  10  to  1  5%  of  a  double  “step 
pe riod  [4,  19]. 

d.  Postural  Control  and  Stability  [22,  23] 

The  foregoing  paragraphs  describe  the  behavior  of  the  lower 
extremities  in  level  walking.  As  far  as  the  upper  part  of  the  body  is 
concerned,  the  center  of  gravity  of  the  trunk  is  imparted  on  an  oscil¬ 
latory  path  as  a  result  of  the  alternate  bipedalism  of  the  gait.  As  one 
leg  serves  as  support  while  the  other  goes  through  swinging  motion, 
there  is  imbalance  of  moments  and  forces  in  the  frontal  and  horizontal 
planes.  These  tend  to  destabilize  the  trunk  from  the  upright  position, 
as  it  resembles  an  inverted  pendulum  mounted  on  a  moving  platform. 
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To  maintain  upright  position  and  minimize  undesirable  oscillations, 
coordination  of  various  muscle  activity  is  needed.  From  a  locomotion 
standpoint,  the  progressive  translation  of  the  body  is  of  main  concern. 
The  associated  unstable  tendency  of  the  trunk  is  an  unfortunate  conse“ 
quence  inherent  in  the  biped  gait. 

e .  A  Hierarchical  Control  System  for  Biped  Locomotion 

From  the  temporal,  phasic  and  stability  properties  of  the  biped 
gait,  it  is  quite  apparent  that  the  locomotion  process  is  involved  indeed. 
It  has  been  postulated  that  the  control  system  possesses  a  hierarchical 
structure,  with  two  or  more  levels  of  decision,  coordination  and 
regulation  for  a  desired  gait  and  adaptation  within  a  given  environment. 
Important  objectives  of  the  control  system  would  be: 

(1)  determination  of  speed,  step  length  parameters  and  direction 
of  progression 

(2)  synchronization  of  the  motion  of  the  two  extremities  to 
achieve  coordinated  phasic  activities  and  temporal  patterns 

(3)  adaptation  to  ground  profile  within  a  double -step 

(4)  maintenance  of  up- right  posture  . 

With  the  exception  of  (i),  the  gait  features  we  have  described 
fit  into  the  objectives  of  the  controller.  From  an  organizational  view¬ 
point,  the  goals  (3)  and  (4)  are  realized  by  a  control  system  under  the 
supervision  of  a  higher  level  controller  connected  to  (1)  and  (2).  A 
schematic  diagram  of  the  locomotion  system  would  appear  as  shown  in 
Fig.  5. 

Such  a  functional  structure  for  the  control  system  possesses  the 
advantages  described  by  Tomovic  [24]: 


SUBCONTROLLERS 


FIG.  5  HIERARCHY  STRUCTURE  FOR  THE  LOCOMOTION  CONTROL  SYSTEM  . 
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(1)  A  minimal  amount  of  supervision  required  on  the  part  of 
controllers  at  the  algorithmic  and  higher  levels  in  executing  well- 
learned  tasks  such  as  normal  level  walking. 

(2)  Flexibility  at  the  dynamic  level  to  adapt  to  a  particular 
environment  within  a  double-step. 

The  present  interest  is  in  finding  out  the  various  actuating 
moments  to  simulate  the  biped  gait.  This  implies  that  the  analysis  is 
performed  at  the  dynamic  level  of  the  hierarchical  system.  The  choice 
of  speed  of  walking,  step-length,  direction  and  other  relevant  parameters 
are  the  goals  of  the  algorithmic  controller,  as  is  the  coordination  of 
muscle  group  activities. 

3.  Energy  Expenditure  and  "Optimality ft  in  Locomotion 

Besides  the  broad  study  of  the  functional  characteristics  of  the 
biped  gait,  the  ergonomics  of  locomotion  forms  another  important  area 
of  research.  An  average  person  walks  over  a  million  steps  in  a  year 
and  certainly  many  times  more  over  the  span  of  his  lifetime.  In  situations 
such  as  level  walking  where  speed  is  not  an  important  factor,  it  is  thus 
very  appealing  to  think  that  the  walking  activity  would  involve  minimiza¬ 
tion  (optimization)  of  some  type  of  energy  performance  criterion.  In 
search  for  such  an  "optimality 11  in  locomotion,  experimental  studies  in 
energy  expenditure,  oxygen  consumption  associated  with  locomotion  have 
been  performed.  Significant  contributions  include  the  works  of  Fenn, 
Elftman  [7],  Ralston  [25],  Cotes  and  Meade  [26],  Earlier  results  were 
obtained  by  Atzler  and  Herbst.  Theoretical  attempts  to  formulate  or 
demonstrate  a  "minimal"  principle  were  made  by  Nubar  and  Contini  [27] 
and  Beckett  and  Chang  [28]. 
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By  far  the  most  interesting  result  has  been  the  measurement  of 
mechanical  work  done  by  the  various  body  members  during  walking  as 
a  function  of  speed,  step-length  and  pace  frequency  parameters.  In 
the  experiments  by  Atzler  and  Herbst  as  reported  in  Elftman  [5],  they 
plotted  energy  expenditure  for  different  speeds  and  step  lengths.  A 
minimum  expenditure  curve  exists  for  the  appropriate  combinations  of 
speed  and  step-length.  In  a  more  recent  study  by  Coates  and  Meade, 
similar  results  were  obtained  and  various  expressions  were  obtained 
from  experimental  data  by  interpolation.  Milner  [20]  and  his  associates 
performed  quantitative  study  of  EMG  activity  of  the  muscles  of  the 
lower  extremities  in  connection  with  the  programmed  electro¬ 
stimulation  problem.  In  particular,  they  showed  that  the  EMG  activity 
is  minimum  when  their  human  subject  is  allowed  to  walk  at  a  pace 
frequency  of  his  own  choice  for  a  given  speed  of  walking.  Whether 
human  locomotion  does  indeed  obey  an  optimality  principle  has  not 
been  proved  by  experiment  or  theory,  but  the  results  gathered  to-date 
do  support  the  suggestion. 

4.  The  Optimal  Programming  Approach 

We  have  briefly  summarized  the  two  important  areas  of  loco¬ 
motion  research.  To  cast  the  biped  gait  in  an  optimization  context, 
the  underlying  hypothesis  is  that  the  locomotion  activity  follows  a 
principle  of  optimality.  In  the  subsequent  analysis,  an  appropriate 
minimizing  performance  criterion  is  developed  which  reflects  or 
measures  the  mechanical  work  done  by  the  muscles  of  the  lower  extrem¬ 
ities  during  walking.  The  system  dynamics  and  auxiliary  constraints 
will  be  developed  on  the  basis  of  gait  information  mentioned  earlier. 
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Numerical  solution  of  the  resulting  optimization  problem  leads  to 
specification  of  the  various  moment  quantities  and  the  simulation  of 
gait  patterns  on  the  basis  of  minimum  mechanical  energy  expenditure. 

This  approach,  using  modern  control  theory,  attempts  to  unite 
the  two  areas  of  experimental  work.  It  differs  from  the  other  studies 
in  that  the  actuating  moments  are  not  predetermined  on  experimental  or 
intuitive  grounds.  The  programming  approach,  if  proved  fruitful, 
opens  a  new  way  to  study  the  locomotion  problem  and  quantify  prosthetic 
design.  To  formulate  the  problem  structure,  specification  is  required 
of  the  following: 

(1)  An  appropriate  mathematical  model  for  simulating  the 
functional  behavior  of  the  locomotor  system. 

(2)  Auxiliary  constraints  -  kinematic  and  dynamic  on  the  basis  of 
gait  information. 

(3)  Initial,  terminal  and  "inflight"  conditions. 

(4)  An  optimality  criterion  which  can  be  extremized  to  yield  the 
actuating  moments  and  other  quantities  of  motion. 
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B .  Mathematical  Modelling  of  the  Biped  Gait 
1  .  Introduction. 

In  a  quantitative  study  of  the  biped  gait,  the  important  moving 
units  are  not  the  various  bones  which  support  the  surrounding  tissues 
but  the  total  mass  of  the  segments  which  rotate  about  their  respective 
joint  axes.  The  central  straight  line  which  extends  between  two  axes  of 
rotation  is  called  a  link.  This  may  be  longer  or  shorter  than  the  parti¬ 
cular  bone  it  represents.  For  example,  the  femur  and  humerus  are 
longer  than  their  respective  links  while  the  tibia  and  radius  are  shorter. 
Functionally,  the  human  body  is  an  articulated  system  of  links. 

Although  the  overall  displacement  of  the  body  is  translatory,  motion  is 
achieved  by  angular  displacements  of  the  links  about  their  respective 
axes.  The  compass  motion  as  mentioned  earlier  is  the  basic  mechanism 
of  biped  locomotion.  The  objective  of  mathematical  modelling  is  to 
describe  the  angular  motions  and  relate  them  to  the  overall  translatory 
process  of  locomotion.  Such  a  mechanical  description  is  separated 
from  the  action  of  the  muscles  and  other  physiological  considerations. 

The  multi -linkage  models  have  been  studied  by  Nubar  and  Contini 

[27] ;  Vukobratovic  and  his  associates  [22,  23,  29];  Beckett  and  Chang 

[28] ;  and  Hill  [30].  Complexity  of  the  models  developed  depends  on  the 
assumptions  made  about  the  gait  and  degrees  of  freedom  allowed  in  the 
motion.  In  a  paper  by  Nubar  and  Contini,  a  system  of  linear  differential 
equations  was  derived  for  small  angle  dynamics.  However,  they  set 
pertinent  derivatives  to  zero  and  studied  only  static  configurations  of  the 
gait.  Marie  and  Vukobratovic  formulated  the  leg  in  level  walking  as  a 
two  link  system.  The  assumptions  on  the  gait  were  that  the  hip  was 
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stationary  and  the  foot  in  continuous  contact  with  the  ground  during  a 
double-step.  There  was  no  mention  of  the  reactional  and  internal  forces. 
In  two  joint  papers  by  Vukob  ratovic ,  Frank  and  Juricic,  the  stability 
problem  of  the  trunk  was  studied  on  the  basis  of  a  compass  gait  for  the 
lower  extremities.  In  the  same  vein,  Beckett  and  Chang  studied  the' 
kinematics  of  swing  phase  motion  under  a  "minimum  energy"  hypothesis, 
but  the  interesting  point  is  that  they  obtained  the  moment  profiles  without 
any  optimization  procedure.  Lastly,  Hill  proposed  a  nine  degrees  of 
freedom  model  for  postural  study.  The  dynamic  equations  are 
apparently  too  lengthy  to  derive.’  Thus  it  appears  that  realistic  efforts 
in  modelling  are  still  lacking.  The  works  done  to-date  share  the 
following  features: 

(a)  They  are  multi -linkage  type  mechanical  models. 

(b)  The  analysis  is  confined  to  motion  in  the  sagittal  plane. 

(c)  Only  the  case  of  level  walking  has  been  treated.  Other 
configurations  such  as  walking  upstairs,  downstairs,  or  over  rough 
terrain  have  not  been  considered  at  all. 

(d)  The  models  have  been  much  simplified  to  render  the  problem 
tractable.  The  result  is  that  they  embody  little  or  none  of  the  described 
properties  of  the  actual  biped  gait. 

2.  Development  of  the  Mathematical  Model. 

The  present  model  emerges  as  a  compromise  between  complex 
high  order  models  and  the  simple  two-link  type.  Under  suitable 
assumptions,  a  high  order  model  is  decoupled  into  two  parts.  The  part 
relevant  to  lower  extremity  is  still  of  a  two-link  type.  However,  appro ~ 
priate  constraints  are  imposed  to  insure  realistic  simulation  of  the  gait 
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patte rns .  Modern  control  theory  can  handle  constraints  effectively 
while  the  classical  methods  cannot.  The  development  of  such  a 
"constrained ”  mechanical  model  takes  place  in  the  following  stages, 
a.  Two  Basic  Configurations. 

Based  on  the  qualitative  features  of  the  biped  gait,  the  actual 
process  can  be  visualized  as  being  made  up  of  two  basic  configurations 
or  modes  of  motion.  Fig.  6  illustrates  the  basic  configurations.  The 
first  is  referred  to  the  double -stand  period  when  both  legs  are  in  contact 
with  the  ground.  One  leg  supports  the  body,  restraining  its  falling 
tendency  as  a  result  of  the  previous  swing  motion.  The  other  leg  starts 
its  deploy  action  through  rotation  of  the  foot  about  its  toes  (more  pre¬ 
cisely,  the  ball  of  the  foot).  We  call  this  the  deploy/restraint  configura¬ 
tion.  The  second  mode  follows  the  first  in  that  the  restraining  leg  now 
takes  up  complete  support  of  body  weight  while  the  other  leg,  previously 
in  deploy,  is  in  its  swing  motion.  We  call  this  the  support/swing  mode. 
Thus,  the  biped  gait  is  an  alternating  sequence  of  these  two  basic 
configurations,  with  the  role  of  the  two  lower  extremities  interchanged 
during  a  double-step. 


Left  leg: 
Right  leg: 


a  double -step 


The  idea  of  two  basic  configurations  introduces  symmetry  into 
the  biped  gait.  It  excludes  singular  behavior  such  as  one  leg  starting  in 
deploy  while  the  other  is  still  in  the  last  stage  of  swing  motion.  In  the 
actual  case,  such  singular  modes,  if  they  exist,  usually  occur  for  a  very 


FIG.  6  TWO  BASIC  CONFIGURATIONS  OF  BIPED  GAIT. 


FIG.  7 


W  =  Weight  of  body 

( R)  =  Restraint 

(S)  =  Support 
( D)  =  Deploy 


[31] 

GROUND  REACTION  PROFILES  (FROM  INMAN). 
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brief  interval  so  that  they  can  be  neglected  in  a  first  study.  In  the 
current  literature,  the  restraint,  support,  and  deploy  portions  are 
collectively  known  as  the  stance  phase  of  activity.  In  our  formulation, 
we  will  refer  to  restraint  and  support  collectively  as  the  "stance11 
portion. 

b.  Derivation  of  the  Model. 

Fig.  8  depicts  a  high  degree  of  freedom  model  for  motion  in  the 
sagittal  plane.  Considering  the  locomotor  system  as  a  whole,  the 
inertial  and  gravitational  effects  of  the  foot  are  small  in  comparison  with 
those  of  the  shank  or  the  thigh.  The  average  weight  of  a  foot  is  about 
1  to  1 .  5%  of  the  body  and  its  dimensions  are  also  small  (approximately 
3  )  as  compared  to  the  thigh  or  the  shank.  Furthermore,  its  behavior  is 
more  or  less  specified  in  a  double-step.  This  is  so  because,  in  the 
stance  position,  the  foot  is  almost  stationary  except  for  a  brief  moment 
following  heel  strike  to  absorb  impact.  In  deploy,  the  ankle  of  the  foot 
describes  a  circular  arc  with  the  center  of  rotation  located  about  the 
ball  of  the  foot.  In  swing,  the  foot  is  locked  withthe  shank  at  an  approxi¬ 
mately  constant  angle.  Such  a  kinematic  behavior  of  the  foot  can  be 
described  by  appropriate  constraints  on  the  thigh  and  shank  variables. 
Dynamically,  the  foot  is  subjected  to  large  ground  reaction  forces  which 
exceed  body  weight  in  the  stance  position.  Characteristic  profiles  of  the 
forces  are  shown  in  Fig.  7.  The  actual  "shape”  of  these  forces  depends 
primarily  on  the  gravitational  and  inertial  effects  of  the  motion  of  the 
upper  part  of  the  body.  The  internal  forces  X,  Y  at  the  ankle  joint  in  turn 
depend  almost  entirely  on  the  ground  reactions.  The  up-shot  of  all  this  is 
that  the  behavior  of  the  foot  can  be  alternatively  described  by  appropriate 


FIG.  8  COORDINATE  SYSTEM  FOR  THE  HUMAN  LOCOMOTOR 
SYSTEM  MODEL. 

NOTATION. :  a,  =  distance  of  center  of  gravity  of  the  thigh  segment  from  hip  joint  (H) 

a2=  . shank  "  "  knee  "  (K) 

1 1  -  length  of  the  thigh  segment 
l2=  .  shank  " 

m1 1  mass  of  the  i-th  segment  ( i  =  1  =>  thigh ;  i  =  2  =>  shank ; 
i  =  t  =>  HAT  section) 

Y,X  ■  vertical  and  horizontal  reactions  at  the  ankle 
Ma=  ankle  moment 

u |  =  effective  moment  for  the  i-th  link 
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constraints  and  equivalent  forces  at  the  ankle.  This  device  reduces  the 
order  of  the  system  dynamics. 

With  the  motion  of  the  foot  considered,  we  are  now  in  a  position 
to  derive  the  mathematical  model  for  the  locomotor  system.  To 
achieve  symmetry  in  the  equations,  let  (h(t),  v(t))  denote  the  horizontal 
and  vertical  position  of  the  hip.  Let  the  variables  <j>  and  y  be  assigned 
to  describe  the  leg  in  the  stance  portion.  The  leg  in  deploy  and  swing 
is  described  by  the  angles  and  x^.  The  head,  upper  extremities  and 
trunk  portions  are  collectively  considered  as  a  single  link  described  by 
the  variable  to.  Note  that  the  hip  position  (h(t),  v(t))  is  a  function  of  the 
variables  cj>,  y.  This  is  also  dependent  on  x^  and  x^  in  the  deploy 
portion,  i.  e  .  : 

h  =  XA  ”  sin  ^  ~  L  +  sin  (Y  "  4  +  L  (l  ) 

v  -  cos  (4  -  40)  ^2  cos  (y  ■  4  40) 


where  (x^,  y^)  denote  the  position  of  the  ankle.  Similar  expressions 
hold  for  the  x^  and  variables  in  the  deploy  portion.  To  carry  out 
the  derivation  by  use  of  Lagrange’s  equations,  expressions  for  the 
kinetic  and  potential  energies  of  the  various  links  are  needed.  The 
centers  of  gravity  for  the  respective  links  are: 

(1)  thigh  in  stance  portion 
=  h  +  Q  1  sin  (4  -  (j>0) 

q>  ■ v " Q  i cos  ‘L 


(2) 


(2)  shank  in  stance  portion 

x  -  h  1  I,  sin  (d>  -  4  )  -  a  0 

y  1  T  To  2 

Y  =  v  -  £,  cos  (4  -  4  )  -  a  0 

7  y  1  1  1  o  2 


sin  (y  "  Q>  +  cj>  ) 
o 

cos  (y  -  4  +  4  ) 


(3) 
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(3)  thigh  in  deploy  and  swing  portions 
x^  =  h  4  a  i  sin  (x^  -  9  ) 

yj  =  v  -  a  j  cos  (xj  -  <j>Q) 

(4)  shank  in  deploy  and  swing  portions 

x~  =  h  4  £,  sin  (x,  -  i  )  -  a  sin  (x~  -  x,  4  4  ) 
Z  L  1  To  Z  Z  1  To 

y  ^  =  v  “  cos  (Xj  "  4o)  “  a  2  cos  (x2  “  xi  +  4q) 


(4) 


(5) 


(5)  trunk  motion  for  all  portions 

/V/ 

x  -  h  h  a  sin  u; 
u  CJ 

y  =  v  4  a  cos  u 
CJ 

The  velocity  of  the  center  of  gravity  for  each  link  is  obtained  through 
differentiation  with  respect  to  time.  From  this,  the  total  kinetic 
energy  of  a  link  is  the  sum  of  translation  component  due  to  motion  of 
C.  G.  and  rotation  about  the  C.  G.  ,  i.  e.  : 


(T)  -  (T  4  T  ) 

1  T-th  link  V  C.G.  Rot.  1  -th  link. 


The  total  kinetic  energy  of  the  link  system  is 

T  =  T  i4  T  4  T  4  T  4  T 
9  Y  xj  x2  u 


(7) 

(8) 


After  manipulation,  rearranging  and  regrouping  terms,  the  total 
kinetic  energy  expression  can  be  concisely  expressed  as 

-  1  A  X  x  1  A  A^  x  1  A  if.  -  X  1  A  ^ 


.  z 


T  =  I  Ao^  +  v  )  +  J  Ai4  +iA2(Y-<i>)  +  k  Alki  +  I  Az^>  "  *1  ) 

•  2  .  f  1  .  !  t 

4  y  A  tJ  4  C,  x,  (ht0  4  vt.  )  4  C0(x0  -  x,  )(-ht A  4  vtA 
2  w  11XZ  1'  Z'Z  1 4  3 


1  *  »  * 

-  -  xj)t6  +  C^ht^  +  vtj)  +  C2(y  -  <j>)(“ht4  +  vt^) 


-  C2<|>(y  “  cj>)t£  +  C^w(h  cos  (J  -  v  sin  u) 


(9) 


21 


where  A  =  2(m,  +  m0)  +  m 
o  1  Z  ( 

A 

i  i  i 

2 


GJ 


1  =  Ij  +  mja  \  +  m2lj 


A2  "  I2  +  m2a  2 


A  -  I  +  m  a 

cj  cj  u;  cj 


C  i  =  m^a  ^  +  m^f^ 


CZ  ~  m2a  2 


C3  m2a  2^1  C2^1 


C  -  =  m  a 

5  cj  u 

sin(cj)  -  cj>o) 

t2  =  cos<4  ■  40) 

4 3  =  sin(  Y  •  4  +  40) 
t4  =  cos ( y  -  4  +  <j>o) 

=  sin(y) 

t&  cos(y) 

Similarly,  the  total  potential  energy  of  the  system  is 


4i  =  sin(x!  -  4>c) 

t 

t0  =  cos(x,  -  <j>  ) 

L,  i  '  o 

1 0  -  sin(x~>  "x,  +  d>  ) 

3  2  1  To 

! 

t  =  cosfx.  “  X,  TO  ) 

4  ^  I  *  o 

! 

t  ^  =  si^x^) 


t6  =  cos(x2) 


g  <AoV  ~  ^lt2  "  (^2t4  "  <^lt2  <^'2t4  +  C^cos  u) 

The  equations  of  motion  are  obtained  by  substituting  the  T  and  V 
expressions  into 


d (£L\. 

dt  \  3q. ) 


9T  ,9V 

9q  3q.  i 


(10) 


(11) 


where  q,  represents  the  angular  variables  and  ML  the  effective  moment 
for  the  appropriate  link.  The  result  is  a  system  of  five  nonlinear, 
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coupled  second  order  differential  equations.  Implicitly,  the  system  of 
equations  can  be  expressed  as 

F  (4> 4? 4; v.  y>  v; *2’ h» v)  =  M  ^ 

4  4 


F(4>  4,  4;  ---  ;h>  v)  = 


(12) 


J 


Fi(4'  4, 4;  ---  ;h,  v)  Mj 

*  •  • 

f2(4, 4- 4;  * hJ v)  -  m2 

V  Fu(u»“>“»h-  v)  =  Mu 

Explicitly,  we  have 
(i)  Leg  in  stance  portion. 

A  (v-  v.  +  h-  h . )  +  A.  cj)  -  A0(  V  -d>)  +  C,  (ht0  +  vt,  )  +  C ,  <i(hf  t0  +  v . t ,  ) 

°  4  4  1  2  1  2  1  1  4 2  4 1 

+  C14Z(-h.t1  +  V  .t2)  -  Cz(-ht4  +  vt3)  +C2(y-  $)(-h.t4  +v.t3) 

4  ’4  4  4 

+  C2(y  -  4)Z(h.t  +  v.t4)  -  C3t6(y  -  2(j>)+C3t5Y(Y“  24> 

4  4 

+  C1x1  (h  .t2  +  V  .tj)  +C1x2(-h.  +  v.t2)  +  C2(x2  -  t4  +  v.  t3) 

cj)  d>  (j>  <j>  <J)  cj) 

2  .  *  .  •  .  v2  - 

+  C  JL  “  x,  )  (h.t0+v.t/1)-f'Cr£>(h/coscj-v.sinu)“C,-£  (h.sincjfv.cosu) 
2  2  1  x  i  3  i  4  5  I  I  7  5  v  ]  1  7 

cp  <J>  cp  cp  cp  cp 


3v 


+  (A°^7+Cltl "  °2t3)g  =  M 


4 


A  (v-  vhh’h  JtAJv-LlCJlh.t.+v.t,)  iCJpEt,  +  v.t0) 
o  y  y  2  1  T  lT  y2  y  1  1 T  yl 

-c2(-ht4  +vt3)  +C2(y  -  4)(-h^t4  +  v^t3)  +C2(y"  4)Z(h^t3  +v^t4) 

V  ;  2  .  t  .  i  2  ■  1  1 

"  C-t,  d) +C -t -4  +  C  t  x\  (h.t0  4-  v  t.) +C,x1(-h.t1  +  v  .t0) 

3  6T  3  5T  1  1  y  2  Y  1  11  y  1  y  2 


t  i 


+  c2(x2-  XjH-ii  t4+v  t3)  +  C2(x2-  x1)2(h  t3+v  t4) 


(13) 
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3v 


+  (AoT7  +  c2t3)g  =  M 


Y 


(14) 


where  h  =  -^-h;  h  =  h:  h,  -  h  and  similar  expressions  hold  for 

dt  dt2  4  34 


v.,  v,  v?  etc.  The  effective  moments  M  j  and  M  are  given  by 


4 


4 

=  u1  +  Mj  +  Y(j?1t1-  l2t  )  +  X(lxt2  +  jf2t4) 


-  ui  T  1Vi4 


M  =  u,  -  M  +  YjLt  -  Xl^t 


Y 


Y 


2  3  2  4 


(15) 

(16) 


Uj  =  Moment  generated  by  muscle  action  about  the  hip  joint 


u2  -  Moment  generated  by  muscle  action  about  the  knee  joint 


M  =  Ankle  moment 
a 


Y  =  Vertical  component  of  the  internal  reaction  force  at  the 
ankle  joint 

X  =  Horizontal  internal  reaction  at  the  ankle. 

Specification  of  the  quantities  X,  Y  and  M^  will  be  considered  later, 
(ii)  Leg  in  deploy  portion. 


!  f 


A  (v*  v  •  +  h*  h .  )  +  A.  x.  -  Ao(x0  -  x,  )  +  C  .  (ht0  +  vt.  ) 
o  x^  x^  it  2  z  I  i  2  I 

+  Clh^x1t2+Y1tl)  +  Cl^l("ix1tl  +  ^x1t2)  '  ^(-ht^+vt^) 
+  C2^x2  "  xl)(“L1t4  +Vi1t3^  +C2^X2  "  xl^  ^x1t3  +  Vx1t4^ 

"  C3l6^2  "  2xl^ +  C3t5x2^x2  ”  2xl^ +Cd^x  tZ  +  Y  Y 


1,+v^.  ,2)+C2(v-^(-h.  <4  +  ^  t  ) 


1 


X1  3 


; ,  2 , 


+  Co(Y  ”  4)  (h‘  t-D  +  v .  t  )  +  Cr  u(h.  cos  W  -  V.  sin  w) 

£  Xj  D  X|  a  *)  X^  X^ 


-  C5u2(h^  sin  u  +  v^^cos  u)  +  (Aq  +  Cjtj  -  C^t^g  -=  M,  (17) 


1 
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Ao(V-  V.  tU-  h  )+A  (^-Xjl+C^jlh  .  tj+i  t|)+Cli|(-h .  tj  +  V^  t2 


•C^tVtjltC^-XjU-i,  t4tv  t'+v  t' 


>  2  *  § 

C,tAx.  +C  t-x.  +C.4(h.  t  +v.  t  )  +  C  d>^(-h.  t.  +v.  t. 
3  b  1  3  3  1  iT  x,  2  x0  1  1T  x0  1  x-,  t 


+  C2(Y-|K-h.^4  +  ^t3)+C2(v-<j)Z(h^t3+v^t4) 

*  <Ao  IT,  +  =  M2  . 

The  effective  moments  and  are  given  by 

Mj  =  u  +  lL  +  -  £zt'3)  +  x'djt^  +  tzt4) 


M  =  u,  -  M  +  Y  Lt,  -  X  Lt. 
2  2  a  2  3  2  4 


(18) 


(19) 


(20) 


'  i  t 

where  ,  Y  ,  X  are  moments  and  internal  reactions  about  the  ankle 
joint  for  the  deploying  leg.  These  quantities  calculated  from  the  tail 

I 

portions  DE  and  D  E  of  the  reaction  profiles  in  Fig.  7. 

(iii)  Leg  in  swing  motion. 

Alh  "  A2^2  "  +Ci(^t2  +vt1)  -  c2(-ht4  +  vt3)  -  C3t6(x2  -  2xj  ) 


+  c3t5^2^2  ”  2i] )  +  (Cjtj  -  C2t3)g  = 


(21) 


A2(x2-  ^1)+Ci(-ht4+vt3)  -C3t6x1  +C3t5i1  +C2t3g  M2 


(22) 


where  M  ^ 


M2  =  u2 


(iv)  Motion  of  the  Trunk. 

A  tj  +  C£-(h  cos  d  “  v  sin  d)  -  Ccg  sin  uj  -  M 
DO  u 


U3) 
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where 

M  =  u  -  u,  (24) 

to  CJ  1  v  ; 

u  ^  =  effective  moment  due  to  muscle  actions  of  the  trunk 
and  upper  extremities. 

Equation  groups  (i)  and  (ii)  form  a  system  of  coupled  equations 
describing  the  motion  in  the  restraint/deploy  configuration.  Note  that 
the  expressions  are  identical  except  for  the  interchange  of  variables  - 
Xj  and  x^  for  cj>  and  y.  This  is  reasonable  as  both  legs  are  in  contact 
with  the  ground  and  the  (h,  v)  -  description  for  the  hip  position  introduces 
symmetry  into  the  expressions.  In  state  variable  language,  the  system 
of  equations  is  equivalent  to  an  eight-dimensional  vector  differential 
equation.  Equation  groups  (i)  and  (iii)  describe  the  support/swing  mode 
of  motion.  The  swing  equations  are  comparatively  simple  and  contain 
no  <{>,  y  variables  explicitly.  Equation  (iv)  describes  the  trunk  motion  in 
the  sagittal  plane. 

c .  Decoupling  and  Simplification  of  the  Model. 

The  foregoing  equations  describe  the  complete  behavior  in  the 
sagittal  plane  and  no  approximations  have  been  introduced.  Motion  of 
the  body  in  the  frontal  and  horizontal  planes  will  introduce  additional 
equations  but  no  modification  in  the  present  ones.  However,  such  a 
Mhigh  order"  model  is  very  complex  to  handle.  Although  the  differential 
equations  can  be  numerically  integrated  as  an  initial  value  problem, 
what  is  at  hand  is  a  multi-point  boundary  value  problem  from  the  optimal 
programming  formulation.  Simplification  of  the  model  is  necessary  to 
render  the  problem  tractable  and  to  obtain  insight  into  the  problem.  The 
major  complexity  of  the  model  is  in  the  coupling  of  motion  between  the 
two  lower  extremities  as  expressed  in  (h,  v)  and  their  various  derivatives. 


The  central  idea  in  simplification  is  to  decouple  the  motion  of  the  two 
legs.  In  level  walking  over  the  range  of  normal  speeds,  pace  frequency 
and  step-length,  the  hip  position  (h,  v)  describes  a  characteristic 
trajectory.  The  horizontal  component,  h,  is  uniform  motion  while  the 
vertical  component  describes  a  double  period  sinusoid  about  some  mean 
value.  Thus,  to  obtain  a  realistic  simulation  of  lower  extremity  motion, 
one  may  introduce  the  auxiliary  constraints 

h  ^  XA  ”  £1  sin^  "  40)  +  *2  sin  ( Y  "  4  +  40)  =  (25) 

v  =  yA  +  *1  cos  (4  “  40)  +  *2 cos  (y  "  4  +  40)  =  g  (t)  (26) 

for  the  leg  in  stance  and  where  f(t)  and  g(t)  are  prescribed  time 
functions  for  the  hip  trajectory.  Similar  expressions  hold  for  the  leg 
in  deploy.  The  auxiliary  constraints  simplify  the  equation  groups  (i), 

(ii)  and  (iii )  somewhat  but  we  still  have  a  high  order  model  as  h„,  h., 

4  Y 

h«  ,  h.  ,  etc.  are  non-zero.  To  obtain  further  simplification,  it  is 
X1  X2 

assumed  that  all  the  partial  derivatives  are  zero,  i.  e.  h*  =  0,  v,  =  0 

4  4 

and  so  on.  Such  a  step  amounts  to  the  consideration  of  the  hip  as  the 
origin  of  a  moving  coordinate  system,  prescribing  the  characteristic 
trajectory  of  the  hip.  Useful  expressions  for  the  hip  position  are 

f(t)  =  vo(t  +  to)  ;  [ft]  (27) 

gU)  =  eo  “  -y sin  ^(t  +  (3  •  T)  ;  [ft]  (28) 

where  v^  =  velocity  of  level  walking 

t  =  constant  based  on  the  initial  configuration. of  the  system 
e^  =  average  height  of  the  hip  above  the  ground 
T  =  period  of  a  double-step 

*3  =  constant  dependent  on  initial  conditions  of  motion. 
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The  device  of  prescribing  a  hip  trajectory  has  been  used  by  Galiana 
[33].  Beckett  [28]  and  Wallach  [34]  in  their  various  studies.  From 
this,  we  obtain  a  decoupled  model  for  the  optimization  problem. 
Simplified  Dynamics. 

(i)  Stance  Portion. 

(Aj  +  2C3t6)4  -  (A^  +  C3t&)y  +  C3t5y(y- 2<j>)  +  (C1t1  +  g) 

-  +  Mq  +  Y(£1t1  -  £zt3)  +  X(£1tz  +  £2t4)  (29) 

-  (A2  +  C3t6)4  +  A2y  +  C3t542  +  C3t6(v  +  g)  =  u2  -  +  Y£2t3 

-  f2t4X  (30) 


( i i )  Deploy  and  Swing  Portions. 

The  form  of  the  equations  is  identical  with  those  of  (i)  except  for 
the  following: 

(a)  The  variables  and  x^  replace  cj>  and  y. 

j  i  i 

(b)  The  terms  ,  Y  ,  X  replace  the  corresponding  terms 
in  (i)  for  the  deploy  portion. 

(c)  In  swing,  all  terms  due  to  ankle  moment  and  reactions 

vanish. 

Because  of  the  similarity  in  the  equations,  we  can  now  adopt  a 
common  notation.  Let  x^  and  x^  stand  for  the  thigh  and  shank  angles 
respectively  (i.  e.  for  (j>  and  y  in  stance;  x^  and  x  ,  in  deploy  and  swing). 
Define  also  x^  -  x^;  x^  =  x^,  then  we  obtain  the  equations  in  first  order 


canonical  form 

x(t)  =  f(x,  u  ;  t) 


(32) 
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whe  re 


d. 


A  ,  ,T 

x  (x1?  x2,  x3,  X^) 


fk  1  =  X3 


x2  -  x4 


<s 


3  =  5  {(R  3  +'ui)  A2  +  t7^Rl  +  u2^ 

^x4  -  ^  {^7(^3  +  ul^  +  t8^R4  +  U2^ 

•  • 

r3  =  ■C3t5x4(x4  "  Zx3>  "  <c1t1  "  C2t3)^  +  g) 


R4  ”  _C3t5X3  ’  C3t6^g  +  g) 

fu^  +  +  Y(£jt^  -  £^t  4-  X(£^t^  +  £^t 4)  in  stance 

^  1  ?  1 

u  i  =  <u^  +  +  Y  (£^tj  ~  £^t^)  +  X  (£^t^  +  £^t^)  in  dePloY 


u. 


*uo 

-  M  +  ’ 

2 

a 

j 

»u. 

-  M*  + 

2 

2 

a 

u2 

V 

6 

_ 

A  ~ 

2  A 
t  -  t.,  = 

2 

8  7 

f 

_ 

A . 

+  2C  t , 

8 

1 

3  6 

t_ 

_ 

A^ 

+*  C  _  t  / 

7 

2 

3  6 

f 

_ 

t 

-  r 

9 

8 

7 

2  3  2  4 

^2^3  ~  ^  ^2^4 


in  swing 

in  stance 
in  deploy 
in  swing 


(33) 


Comments . 


Using  the  dynamics  in  first  order  form,  one  can  study  the  motion 
of  the  two  lower  extremities  by  considering  only  the  sequential  behavior 
of  a  single  leg.  Although  the  expressions  are  similar  to  those  obtained 


29 


in  two-link  modeib,  we  have  placed  into  perspective  that  the  simplified 
model  is  obtained  through  decoupling  of  a  high  order  one.  In  the  model 
by  Beckett  and  Chang,  only  the  kinematic  aspects  of  the  motion  are 
considered  by  neglecting  all  reactions  and  ankle  moment  contributions. 
Wallach's  [34]  paper  on  knee  prosthetics  also  neglected  the  forces  and 
moment.  In  the  paper  by  Galiana  [33],  the  ankle  moment  term  was 
missing  as  a  result  of  the  !,fixed~foot  ”  assumption.  The  decoupled 
equations  are  identical  to  those  derived  by  Bresler  and  Frankel  [13] 
using  the  D'Alembert’s  principle  of  dynamic  equilibrium  for  their 
experimental  study.  Their  important  results  indicate  that  the  ground 
reaction  forces  in  the  stance  portion  are  the  most  significant  quantities 
and  gravitational  and  inertial  effects  of  the  thigh  and  shank  variables 
are  small  compared  to  these.  Thus  it  appears  that  the  simplified  model 
should  yield  realistic  results  provided  appropriate  consideration  of  the 
foot  motion  is  taken. 

e .  Kinematic  Constraints. 

In  the  deploy  portion^  the  ankle  of  the  foot  describes  a  circular 
arc  about  the  ball  of  the  foot.  With  the  hip  motion  prescribed,  the 
angles  and  are  no  longer  independent  but  subjected  to  an  equality 
constraint  relationship.  Suppose  that  deploy  starts  at  time  t^  and  ends 
at  t^.  Vertical  displacement  is  summed  to  yield 

eQ  ='g(tl)  +  ^  j  cos(xj  (tj )  -  4q)  +  £2cos(x2(t1)  -  x^tj)  +  <j>Q) 

+  dsina(t,  ) 

as  shown  in  Fig.  9.  The  angles  x^  (t  ^  ,  x^(t  ^  )  and  a  (t^  )  are  specified. 
Similarly,  for  <  t  ^  t^,  we  have 

eQ  =  g(t)  +  i j  cos(xj  (t)  +  ^cos(x^(t)  -  x^  (t)  +  <|>o)  +  dcosa(t)  . 


(Deploy)  (Swing) 

FIG.  9  KINEMATIC  CONSTRAINTS  IN  DEPLOY  AND 
SWING  PORTIONS. 


FIG.  10  KINEMATIC  CONSTRAINTS  IN  STANCE  PORTION. 
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Eliminating  between  the  two  expressions,  we  have 

d  sin  a  (t)  =  ej  -  g  (t)  -  =  S1  =  yA  (34) 

el  =  g  (tj)  +  (ijt2  +  i2t4)-  +  d  sin  a  (t^). 

Summation  of  horizontal  displacements  gives 

vQ(ti  +  tQ)  +  i1sin(x1(t1)  -  <j»o)  -  ^2sin(x2(t1)  -  x^)  +  (j>Q) 

+  d  cos  a  (tj)  =  vQ(t  +  to)  +  £1sin(x1(t)  -  <j»Q)  "  *2sin(x2(t)  -  x^t)  +  cj>o) 

+  d  cos  a  (t) 

or  d  cos  a  (t)  =  -  VQt  -  4-  =  S2  =  xA  (35) 

e2  =  v  4-  (f^t^  -  ^2^3^  +  ^  COS  a  ^1^ 

Eliminating  a  (t)  between  and  S^,  we  obtain  the  desired  relationship 

Sd(x;t)  =  +  ^2  “  d2  +  (ej  *"  1>)Z  +  (e2  -  v^t)Z  +  2^4^ 

-2(ei  -  I)(f1t2  4-  l2t4)  -  2(e2  -  vot)(f1t1  -  l2t3)  =  0.  (36) 

At  time  t2?  the  angle  6  between  the  shank  and  foot  approaches 
a  limiting  maximum  value,  6  max.  During  swing  motion,  6=6  max 
(i.e.  foot -locking)  while  the  toes  of  the  foot  are  kept  above  the  ground. 

From  Fig.  9,  we  thus  have  the  inequality  constraint 

Ss(x;t)  =  g(t)  4-  £^t2  *  ^2^4  +  dcos|3(t)-e^^0  (37) 

(3  =  77  -  6  max  -  (j>0  +  Xj(t)  -  x2(t)  (38) 

for  t2  ^  t  ^  t^,  being  the  end  of  swing  (i.  e.  -  T  4-  t^  ; 

T  =  period  of  a  double- step) . 

In  addition  to  state  variable  inequality  constraint  (37),  the  angle 
x2(t)  has  to  be  positive  (x2  >  0)  during  the  swing.  However,  the  multi¬ 
linkage  mechanical  model  does  not  take  into  account  such  a  restriction. 
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This  aspect  of  model  limitation  does  not  appear  to  have  been  noticed  in 
previous  investigations.  A  useful  device  to  alleviate  the  difficulty  is  to 
adjoin  a  "soft"  penalty  [35]  to  the  performance  criterion  under  optimi¬ 
zation"^,  i.  e.  define 


J 

s 


dt 

x2(t) 


(39) 


for  a  sequence  of  yq  >  >  Y ^ - >  Yg— ’ ►0^-  The  penalty  becomes 

effective  as  x^— *  0_^_  but  smaller  as  x2(t)  becomes  more  positive.  A 
restriction  is  that  in  iterative  computation,  the  trajectory  x^(t)  has  to 
remain  feasible,  i.  e.  x2(t)  >  0  Vt  e  [f?,t^].  Intuitively,  the  factor 


- jr,  is  equivalent  to  a  nonlinear  spring  placed  about  the  knee  joint  and 

the  parameter  y^  measures  its  stiffness. 

In  the  stance  portion,  the  foot  is  more  or  less  stationary  except 
for  the  brief  moment  following  heel  strike.  The  dynamic  behavior  is 
difficult  to  describe  because  a  rigid  link  cannot  describe  the  behavior  of 
a  "round"  foot.  Fig.  10  illustrates  the  roll-over  situation  where  the 


weight  bearing  leg  rotates  about  its  ankle  (0  )  -  the  primary  rotation 
center.  The  pressure  center  where  the  ground  reaction  forces  are 
acting  is  gradually  transferred  from  the  heel  towards  the  toe.  The 
effective  kinematic  constraints  in  this  period  are  that  the  ankle  (x^,  y^) 
coordinates  follow  a  prescribed  trajectory,  (q(t),p(t)).  In  a  rough 
analysis,  q(t)  and  p(t)  can  be  taken  as  constant.  Thus,  we  have 


S[s(x;t)  =l(t)  +  Ijt 2  +  <2 14  -  eo  +  p(t)  =  0  (40) 

S2S(x;t)  =  VQ(t  +  tQ)  +  l1t1  -  ^3  ‘  q(t)  =  0  (41) 

-}  Jq  has  yet  to  be  defined. 
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The  particular  choice  of  p  and  q  profiles  has  to  depend  on  experimental 
data.  Fig.  11  shows  the  typical  profiles.  Note  that  they  are  almost 
constant  except  for  the  initial  restraining  moment. 

f.  Quasi -Static  Estimation  of  Ankle  Moment  and  Reactions. 

Having  developed  the  dynamic  equations  and  kinematic  constraints, 
the  model  is  incomplete  until  the  ankle  moments  and  reaction  forces 
are  known.  Fig.  12  shows  a  free  body  diagram  of  the  foot.  The  forces 
f^  and  f  represent  the  vertical  and  horizontal  components  of  the  ground 
reactions  due  to  weight  bearing  and  trunk  motion.  These  forces  have 
been  measured  by  various  investigators  using  force  plates  and  Fig.  7 
shows  the  typical,  normalized  profiles  over  a  wide  range  of  level  walking 
speeds  and  step-lengths.  To  describe  the  foot  motion,  we  write  the 
equations  for  translation  motion  for  the  center  of  gravity  (O^)  and 
rotation  about  O^.  On  summing  the  vertical  forces,  we  have 


or 


Y  =  fy  -  MFg  -  (MFyF) 


(42) 


Similarly,  summing  the  horizontal  forces  gives 


M  ’x^  =  f ,  “  X 
F  F  h 


or 


X  =  f  -  M  x 
h  F  F 


(43) 


Taking  moments  about  and  rearranging  terms  yields 


-  (MFyF)(yF  -  yA)  +  Ij-d 

^  -  2  2 

1  =  Moment  of  inertia  of  a  foot  about  its  C.G.  (~  10  slug-ft  ) 

F 


(44) 


I  l 

q(t) 


restraint 


support 


i+d  cos(t|) 

J  =  step  length 


FIG.  1 1  TYPICAL  PROFILES  FOR  ANKLE  COORDINATES  IN  STANCE 
PORTION.  (FROM  BRESLER  AND  FRANKEL) . 


Op  *  pressure  center  ( xp ,  yp ) 

Op  -  center  of  gravity  of  the  foot 

CxFtyF) 

°a  =  onkle  {h>h] 

X,Y  *  reactions  (internal)  (a)  the  onkle 

fv,fh  *  vertical  and  horizontal 

components  of  ground  forces 


FIG.  12 


FREE  BODY  DIAGRAM  OF  THE  FOOT. 
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It  is  obvious  that  the  ankle  reactions  and  moment  depend  on  the  ground 
reactions.  Since  fy  »  (Mpg);  fy»  (MFYp)»  fh  »  (MF^p);  and  (1-p.a  ) 
is  small  because  Ip  is  small,  we  can  obtain  a  simple  "quasi  -  static  M 
estimate  of  the  dynamic  quantities  by  neglecting  the  inertial  and  gravita¬ 
tional  terms,  i.  e.  : 


Y  ~  f 

v 

x~fh 

Ma  -  £v(xp  '  XA>  +  fh(yA  '  V 

The  ankle  coordinates  (x^,  y^)  are  specified  through  the  kinematic 

constraints  derived  in  the  preceding  section.  As  far  as  the  pressure 

center  (x^,  Y  )  is  concerned,  we  have  y^  =  0  throughout  the  stance 

portion.  In  the  report  by  the  Ebe  rha  rt -Inman  group  [12],  force  plate 

data  reveals  an  approximately  linear  transfer  characteristic  in  the 

stance  portion,  i.  e.  : 

x  ~  a  t  +  b 
P 


(45) 


(46) 


where  a  ,  b  are  parameters  dependent  on  walking  speed,  foot  length 
and  other  parameters.  During  deploy,  x^  is  fixed  at  the  ball  of  the 
foot.  With  the  dynamic  quantities  specified,  the  multi  -  linkage  model 
with  auxiliary  constraints  is  now  complete. 

3.  Discus  sions . 

We  have  derived,  through  a  succession  of  steps,  an  appropriate 
model  for  the  optimization  problem.  It  is  a  two-link  model  with  kine¬ 
matic  constraints  and  ankle  forces  and  moments.  Fig.  13  is  a  flow 
chart  illustrating  the  development.  Note  that  models  at  different  levels 
of  detail  can  be  obtained  by  allowing  additional  degrees  of  freedom  and 


FIG.  13  FLOW  CHART  OF  MODEL  DEVELOPMENT. 
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complexity. 

An  inherent  drawback  of  the  high  order  model  is  that  the  dynamic 
variables  chosen  for  description  are  always  coupled.  This  implies  the 
matrix  inversion  is  necessary  to  convert  the  equations  to  first  order 
form.  For  models  with  more  than  two  degrees  of  freedom,  matrix 
inversion  with  analytic  expressions  soon  becomes  too  messy  to  handle. 

The  use  of  symbolic  manipulation  language  such  as  FORMAC  appears 
indispensible  in  the  study  of  complex  models  in  locomotion  and  bio¬ 
engineering  in  general. 

The  use  of  a  prescribed  hip  trajectory  and  kinematic  constraints 
for  ankle  motion  is  a  useful  device.  Locomotion  under  different  conditions 
such  as  up  an  inclined  surface  or  over  other  terrain  only  changes  the 
form  of  the  prescribed  functions  and  not  the  other  expressions.  Since 
modern  control  theory  can  handle  readily  auxiliary  constraints,  it 
appears  very  promising  in  the  study  of  locomotion  under  various  environ¬ 


mental  conditions. 
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C .  On  an  Energy  Performance  Criterion  in  Human  Locomotion 

1.  Introduction. 

Having  developed  an  appropriate  model  for  the  optimal  programming 
problem,  the  next  step  is  to  derive  a  suitable  performance  criterion.  As 
mentioned  earlier,  Atzler  and  Herbst  [5]  and  Cotes  and  Meade  [26] 
independently  obtained  experimental  results  which  show  a  minimum 
expenditure  of  mechanical  work  for  appropriate  combination  of  speed, 
step-length  and  pace  frequency.  Theoretically,  Nubar  and  Contini  [27] 
were  probably  the  first  and  only  thus  far  to  propose  a  minimal  principle 
for  human  locomotion. 

They  hypothesized  that  energy  expenditure  under  several  simul¬ 
taneous  forms  (mechanical,  heat  and  chemical)  is  associated  with 
muscular  activity.  The  muscular  effort,  an  aggregate  quantity,  is 
defined  quantitatively  in  terms  of  the  effective  moments  which  actuate 
the  locomotor  system.  The  locomotion  process  is  thought  to  occur  in 
such  a  manner  as  to  minimize  some  "effort*1  function  E^  (performance 
criterion)  consistent  with  imposed  constraints  -  physiological,  physical 
and  geometrical.  The  effort  function  they  proposed  is  of  a  quadratic 
variety,  i.  e. 

V  A  2 

Er  =  /  C.  u  .  (t)  At 

f  L,  }  } 
j 

where  u  .(t)  represents  the  effective  moment  acting  on  the  j-th  link; 

A  J 

Ck  being  an  effectiveness  weighing  factor;  and  At  is  the  time  interval 
over  which  optimization  is  performed. 

Their  work,  however,  has  two  important  shortcomings.  At  that 
time,  optimal  control  theory  and  numerical  computation  were  not 
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developed  extensively  and  the  proposed  problem  was  not  seriously 
tackled.  A  second  and  more  fundamental  question  is  that  while  the 
effort  function  being  chosen  is  selected  on  the  basis  of  mathematical 
tractability  and  physical  appeal,  there  is  no  apparent  connection  to  the 
physiology  of  the  muscle  actuating  system.  Our  present  objective  is 
to  derive  a  simple  but  physiologically  based  performance  index. 

2 .  Development  of  the  Performance  Criterion. 

Mechanical  energy  expenditure  is  the  quantity  of  interest  for  the 
derivation.  In  general,  mechanical  work  done  (by  or  on)  can  be 
evaluated  in  two  ways.  One  measure  is  the  time  integral  of  the  instan¬ 
taneous  power  which  is  the  product  of  the  moment  about  a  joint  and  the 
angular  velocity  of  the  limb  with  respect  to  that  joint.  The  other 
measure  is  the  integral  of  the  product  the  force  generated  by  the 
muscle  in  causing  rotation  and  its  velocity  of  shortening.  The  latter 
method  is  used  as  it  involves  naturally  the  muscle  mechanics. 
Specifically,  the  external  characteristics  of  the  muscle  are  used, 
a .  Assumptions . 

The  neural- viscoelastic  model  for  the  muscle  is  based  on  the 
following  set  of  experiments  which  have  been  performed  on  both 
isolated  muscle  fibers  and  the  human  muscle  in  vivo 

(i)  tension  (force)  -  velocity  curves  by  A.  V.  Hill  [36,  37] 

(ii)  length  -  tension  curves  [36,  37] 

(iii)  EMG  characteristics  by  Bigland  and  Lippold  [46] 

Fig.  14  illustrates  these  properties.  From  these,  the  force  generated 
by  a  muscle  is  a  function  of  neural  stimulation  Z,  muscle  length  t  and 
velocity  v.  The  functional  form  is  given  as  [38,  39] 


MUSCLE  OUTPUT  FORCE 


SHORTENING  VELOCITY 


TENSION 


[36,37] 

FIG.  14  EXPERIMENTAL  CHARACTERISTICS  OF  SKELETAL  MUSCLES 


FIG.  15  AGONIST  /  ANTAGONIST  MUSCLE  SYSTEM. 


where  Fq(!)  maximum  isometric  force  at  length  i 

Z  =  neural  stimulation  level  (0  ^  Z  ^  1) 

v  maximum  velocity  for  isotonic  case  (i.  e.  when  F  0) 

m 

b  -  constant  (typical  values  between  0.  25  and  0.  3). 

Thus,  F  is  a  function  of  three  variables,  f ,  v,  Z.  By  holding  one 
variable  constant,  the  muscle  length  H  in  this  case,  F  can  be  graphically 
portrayed  as  a  surface.  This  we  call  the  characteristic  surface  of  the 
muscle . 

The  derivation  is  confined  to  one  joint  agonist/antagonist  muscle 
system.  The  actual  muscle  system  involved  in  almost  any  complex 
limb  motion  is  seldom  that  simple.  We  assume  that  the  same  principles 
hold  for  each  agonist/antagonist  pair  involved,  and  thus  the  discussion 
represents  an  average  behavior  of  all  pairs  contributing  to  the  actual 
limb  motion.  The  human  locomotor  system  is  equipped  with  an  abun¬ 
dance  of  two  joint  muscles;  however,  Steindler  [4],  Eberhart  et  al  [10] 
have  shown  how  two  joint  muscles  can  be  functionally  decomposed  into 
the  one  joint  variety, 
b.  Derivation. 

Consider  a  limb  being  acted  upon  by  an  agonist/antagonist 

muscle  system  (Fig.  15).  With  the  two  moment  arms  approximately 

equal,  there  is  no  net  rotation  of  the  limb  when  each  muscle  exerts  an 

equal  force  F  =  Z  F  (I  ).  The  presence  of  static  forces  is  to  maintain 
^  o  o  o  o 

joint  stability  and  body  posture.  When  one  muscle  exerts  more  force 
than  the  other,  a  net  rotation  results.  The  amount  of  mechanical  work 
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done  by  the  muscle  actuating  system  is  given  by 
W 


=  \  (F  t  6  F  )  •  v  dt  +  I  (F  +  6  FJ  •  v  dt 
J  o  s'  s  J  o  r  e 


(2) 


1  A 


F  (v  +  v  )dt  +  \  6  F  •  v  dt  +  \  6  F.  v  dt 
o'  s  e  J  s  s  .  £  e 


where  6  F  and  0  F  are  the  incremental  forces  above  the  static 
s  e 

A 

component  F  in  shortening  and  lengthening  respectively.  The  corre¬ 


sponding  shortening  and  lengthening  velocities  are  v  (>  0)  and  v  (<  0). 

s  e 


In  normal  locomotion,  it  is  taken  that  jv 


v  !  and  so  the  static 


dependent  term  F^(v^  +  v^)  ~  0,  contributing  little  or  nothing  to  W. 
Muscle  in  Shortening. 

From  (1),  the  "shortening 11  force  is 

1  -  JjL 

V 

F  =  ZF  (£)  - — - 

S  O  V 


(3) 


1  - 


m 


For  normal  locomotion  activity,  F^  is  not  significantly  different  from 
the  static  component  so  that  it  can  be  satisfactorily  represented  as 


5  F  =  F  -  F 
s  s  o 


and 


Z  -  Z  )  + 
o 

o 


(l  - 


o 


( small) 


The  fact  that  the  linear  terms  constitute  the  main  contribution  is 
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supported  by  experimental  results  obtained  in  normal  walking. 

(i)  The  stimulation  level  Z  appears  linearly  in  (1)  and  so  there 
is  no  approximation  in  (4). 


(n)  The  term  due  to  length 


small  compared  to  the  other  two;  consequently  it  is  omitted  Results 
by  Elftman  on  the  soleus  muscle  indicated  that  the  muscle  length 
change  is  small  during  the  shortening  cycle  (Fig.  16).  Also,  from  the 
length-tension  diagram  the  shape  of  the  curve  is  almost  flat  in  the 
vicinity  of  the  rest  length.  In  the  EMG  experiment  by  Galiana  [33], 
satisfactory  results  were  obtained  under  a  length  independence  assumption 
(iii)  The  effect  of  shortening  velocity  v^  on  F^  is  essentially 
linear  in  the  range  of  interest.  The  expansion  generates  a  linear 
approximation  to  the  hyperbolic  curve  of  Hill.  It  is  interesting  to 
observe  that  the  stimulation  level  Z  is  multiplicative  rather  than  auditive, 
so  that  it  controls  the  damping  as  well  as  the  isometric  force  effectively. 
Physically,  the  linear  approximation  replaces  a  velocity  dependent 
damper  by  a  constant  one  whose  magnitude  depends  on  the  isometric 
conditions.  Stark  [41]  and  Houk  [42]  had  used  linearization  in  their 
modelling  of  muscle  actuators;  their  results  appear  satisfactory. 

By  neglecting  the  length  contribution,  we  have 
fiF  =  F  -  F 


s 


s 


o 


s 


where  B  =  Z  F  (£  )(1  Ir)  /  v 

m  o  o  o  b  m 


The  next  crucial  step  is  to  obtain  a  relationship  between  stimu¬ 


lation  increment  A Z  =  Z  -  Z^  and  velocity  v  .  Bigland  and  Lippold 


s 


RESULTS  ON  THE  SOLEUS  MUSCLE 


5  cm 
4  cm 


FIG.  16  VELOCITY-LENGTH  CHARACTERISTICS  (from  ElftmanC7:i) . 


FIG.  17  CHARACTERISTIC  SURFACE  OF  A  SKELETAL  MUSCLE. 
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showed  that  for  approximately  constant  and  AZ  are  linearly 


related  for  small  velocities.  This  is  easy  to  verify  from  (3): 
F  ^ 

(1  +  b)  — —  cC  v  fori  v 


AZ 


F¥(1  +b>7S-<»c  v8  for'— ^ 
o  m 


i«  1 


(6) 


m  / 


In  the  normal  range  of  locomotion  activity,  a  linear  relationship 

between  AZ  and  v^  still  holds  but  with  a  slightly  larger  slope  dependent 

on  the  maximum  value  of  F  .  The  "equivalent  linearization"  is  clear 

from  the  characteristic  surface  as  illustrated  in  Fig.  17.  The  surface 

is  a  graphical  representation  of  the  functional  relationship  (1)  for  f  - 

constant.  As  the  muscle  length  changes  are  fairly  small  in  normal 

locomotion,  the  locus  of  F  lies  within  two  characteristic  surfaces 
5  s 

lying  close  together.  For  F^  =  constant,  increasing  stimulation  AZ  lies 
on  paths  AC  or  DB  on  the  characteristic  surface.  With  F  varying 

A  A  ^  A 

from  F  to  (F  )max  ~  1.  IF  (estimated  from  Elftman’s  result),  the 
o  s  o 

r\ 

locus  is  the  curve  AB  on  the  surface.  When  projected  onto  the  Z  -  v 
plane,  we  have  an  almost  linear  segment  A’B’.  Thus,  we  have  an 
operating  relationship  of  the  form 


A  A 

AZ  =  k  ((F  )max:  F  )  *  v 
s  s  ’  o  s 

Note  that  if  one  uses  the  line  A’C!  instead  (i.e.  no  correction  for  F 


(7) 


changing),  the  error  is  approximately  10%  for  the  range  of  F  variation. 


Combining  (5)  and  (7),  we  have 

6  F  =  (F  k  -  B  )  •  v 
sos  m  s 

s 

Since  (F  -  F  )  >  0  for  a  shortening  muscle,  then  (F  k  -  B  )  >  0. 
'so  ’  o  s  m 

s 


(8) 


The  work  done  by  a  muscle  in  shortening  is  then 


4  7 


W 

s 


i‘v.  -  * 

(F  k  -  B 
Jos  m 

t 

o 


) 

s 


dt 


u  (t)  dt 
s 


(F  k  -B  )  •  d2(t) 
os  m  s 

s 


2 


r  (t)u  (t)dt 

s  s 


(9) 


where  u  (t)  is  the  moment  generated  by  the  shortening  muscle  and 


d  (t)  the  moment  arm. 
s 


2 


rs(t) 


A 


1 


(F  k  -B  )d2 
os  m  s 
s 


(10) 


Muscle  in  Lengthening. 

When  the  other  member  of  the  agonist/antagonist  pair  is  being 

stretched,  mechanical  work  is  done  on  it.  The  Hill’s  force-velocity 

curve  holds  for  the  lengthening  muscle  up  to  a  maximum  allowable 

force  which  is  roughly  twice  F  (£  ).  For  larger  forces,  the  Golgi 

tendon  organs  reflexly  cause  the  muscle  to  yield. 

Following  the  same  procedure  as  in  the  preceding  case,  the 

force  in  lengthening  is  given  by 

6  F  =  F  -  F 

i  a  o 

=  F  AZ  -  B  v  (11) 

o  v  7 

B 

where  v-  <  0  and  _ ^£  ~  ,  according  to  Stark,  Galiana.  Thus, 

B  ^ 

m 

s 


the  linear  approximation  to  the  characteristic  surface  is  a  piecewise 
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one  with  the  damping  during  lengthening  (B  )  much  larger  than  that  in 

m£ 

shortening.  Fig.  18  illustrates  that  below  saturation,  AZ  and  are 
linearly  related,  i.  e 


AZ  =  k,((FJ  :  F  )v„ 
£  £  max*  o  £ 


The  work  done  on  the  muscle  is 


W  =  - 

e 


A  2 
(F£  -  Fo> 


(B  -  F  k„) 
m£  °  1 


■dt 


(12) 


(2) 


u^(t)dt 


(B 


A 

-  F  k.) 

mi  °  £ 


df(t) 


r£(t)u£(t)dt 


(13) 


where  r^(t)  = 


- 2 —  and  d^(t )  =  moment  arm  for  the 


lengthening  muscle  d  (t)). 

Total  Mechanical  Work. 

The  total  work  done  is  the  sum  of  lengthening  and  shortening 

contributions 

W  =  +  W 

£  s 

ff 

=  t  1  (r  u2  -  r  u^dt. 

2  J  v  s  s  £  £ 

t 

o 


(14) 
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The  net  moment  acting  on  the  limbs  is  given  by 


u 


F 

s 


-  F 

d 


s 


£ 


A  A 

(F  -  F  )  -  (F  -  F  ) 

S  O  6  O 

d 

s 


Now, 

But, 


~  us  "  U*  for  tf'*1 

s 

/  i2  2  .  2  , 

=  (US  V  1  US  +U*  '  insUl 


U  , 


U  +  U, 


2<US+V 


(15) 


u 


W 


k  ,  2  Z 

=  y  (u  -  uJ  ;  k  =  constant  factor 


2  v  s 

ij  (rsuf  'r(uj)dt 

t 


-  constant 


u^(t)dt 


(16) 


t 

o 

Hence,  in  the  normal  range  of  activity,  the  sum  total  of  mechanical 
energy  expenditure  by  the  muscle  activating  system  is  proportional 
to  the  integral  of  the  square  of  the  net  moment. 

3.  JDi  scus  s  ions . 

In  the  present  derivation,  mechanical  energy  expenditure  by  the 
muscle  system  is  the  physical  quantity  being  identified.  Although  other 
quantities  such  as  total  energy  expenditure  by  the  muscles  or  total 
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"active  states  11  of  the  muscles  have  been  suggested  as  bases  for  a 
performance  index,  it  appears  not  straight-forward  to  develop  these. 
The  mechanical  work  done  by  (and  on)  the  muscles  has  been  of  most 
concern  in  experimental  work  done  to-date. 

As  to  the  form  of  the  effort  function,  generalization  can  be 
made  of  the  form 


o 


where  p  is  some  appropriate  integer;  u.  the  net  moment  acting  on  the 
i-th  link  and  the  weighing  factor.  Physiologically,  justification  for 
one  particular  choice  may  be  more  difficult  than  for  another.  For 
the  case  p  =  0  and  t^  treated  as  unspecified,  the  criterion  leads  to  a 
"minimum  time"  problem  and  appears  applicable  to  fast  walking.  The 
p  -  1  case  corresponds  to  the  "minimum  fuel"  problem  in  optimal 
control  theory.  For  p  =  2,  we  have  the  quadratic  criterion  considered 
in  the  preceding  sections. 
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D.  Optimization  in  Biped  Locomotion 


1.  Problem  Formulation  [431 

The  constrained  mechanical  model  and  performance  criterion 
form  the  core  of  the  optimization  problem  for  biped  locomotion. 
Based  on  the  simplified  model,  simultaneous  motion  of  the  two  legs 
is  reduced  to  the  sequential  behavior  of  a  single  leg.  Let  deploy, 
swing  and  stance  portions  describe  in  order  the  phasic  activities  in  a 
walking  cycle  (i.  e.  double  -  step) .  Mathematically,  the  optimization 
problem  can  be  stated  as  follows: 

Let  t^  be  the  time  required  for  a  double-step  and  t^  and 
partition  [0,  t^]  into  the  respective  phasic  portions.  It  is  proposed 
that  level  locomotion  is  realized  by  programming  the  hip  and  knee 
moments  (u^  u^)  so  as  to  minimize  the  quadratic  criterion. 


ft 


1 


(  r  |  u  j-fr^i^Jdt  - 


\urTz  2 


J  0 


1 


z 


for  rp  r^  >  0  ;  0  <  <  t^  <  t^.  and  subject  to 

a.  dynamic  equations. 
x(t)  =  f(x,  u;  t) 


(2) 


i.e.  ,  equations  (32),  (33)  of  the  modelling  section. 


b.  kinematic  constraints. 


for  0  ^  t  ^  t 


1 


(3) 
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S‘S(x;  t)  -  g(t)  +  +  ^2l4  +  d  cos  P  '  e0 


P  =  37-6  "  <j>  +  x,  - 

max  *o  1 


X,- 


for  t^<t^ 


(4) 


x2  >  ° 


S^°  (x;t)  =  g(t)  +  i1t2  +  f,t4  +  p(t)  -  eQ  =  0 


r  for 


Srzs  (x;t)  =  vo(t  +  to)  +  tj  -  £2t3  -  q(t)  =  0, 

c •  initial  and  terminal  conditions 
x(0)  =  x^  (specified) 

^X(x( ^[Sd(x;t)]t=t  -  (Sj  +  S2  -  d\  =0 

y2(x(tl);tl)  =  x2(tj)  -  xjdj)  +  4o  +  6  max  -  Q{1 

These  are  the  terminal  hypersurface  constraints  for  deploy. 
x(t2)  =  X2  “  d  (specified) 

x(t^)  =  x(0)  (repeatability  in  walking) 


(5) 

(6) 

(7) 

(«) 

(9) 


Note  1  . 

Conditions  (6),  (7)  and  (8)  describe  the  configuration  of  the  leg 
at  heel  off,  toe  off  and  heel  strike  respectively.  In  our  study,  the 
values  x(0)  and  >l( t ^ )  are  specified  from  experimental  data,  as  they 
together  depict  the  deploy/re straint  configuration  at  the  starting 
instant.  Such  a  specification  appears  more  natural  than  Becketths 
case  where  x(t^)  remains  unspecified. 

Note  2. 

A  hypersurface  type  of  constraint  is  used  to  describe  the  toe  off 
condition  at  the  end  of  deploy.  Although  the  state  equality  constraint 
S  (x;  t)  is  binding  throughout  the  deploy  portion,  the  additional  condition 
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Y  ^  =  0  is  a  device  used  to  eliminate  possible  jumps  for  the  influence 
equations  in  numerical  simulation.  The  ^  ^  =  0  condition  describes 
"foot-locking "  for  the  subsequent  swing*  and  6  max,  a(t^)  are  set  at 
110C  and  70°  for  simulation  [10]. 

Note  3. 

The  condition  (9)  describes  repeatability  of  motion  in  normal 
level  walking. 

Note  4. 

From  the  kinematic  constraints  (3)  describing  deploy  motion  and 
Fig.  12,  we  can  write 

=  el  -  g(t)  -  tl  t2  -  lz  t4  =  d  sin  a 

S  =  e0  -  v  t  -  £  t  +  £  t_  =  dcosa 
2  2  o  11  2  3 

where  d  =  effective  length  of  the  foot  between  the  ankle  joint  and 
the  ball  of  the  foot  (i.  e.  center  of  rotation) 


a  (t) 


d  (t)  = 


S1S2 


2.  Characteristics  of  the  Problem  Formulation 

a.  It  is  a  programming  approach  as  opposed  to  the  direct  approach 
taken  in  other  studies. 

b.  The  formulation  results  in  a  multi-arc  programming  problem 
with  three  stages  on  account  of  the  varied  nature  of  the  ankle  reactions 
and  moments;  kinematic  constraints;  initial  and  terminal  conditions. 


56 


c.  Computationally,  the  multi-arc  problem  is  decomposable  into 
three  problems  for  solution.  Continuity  in  the  state  variable  histories 
(angular  displacements  and  velocities)  is  ensured  by  matching  the 
trajectories  between  two  adjacent  stages. 

d.  The  programming  formulation  is  very  flexible.  Although 
emphasis  has  been  on  the  human  case,  the  formalism  and  approach 
are  certainly  applicable  to  practical  designs  such  as  above -knee 
prosthetics  and  remotely  controlled  manipulators  utilized  in  a  hostile 
environment . 

e.  The  splitting  up  of  the  multi-arc  programming  problem  for 
independent  numerical  solutions  links  to  the  important  concept  of 
sub -optimality  in  our  solution.  For  the  present  simulation,  the  inter¬ 
mediate  times  t-^  and  t^  are  specified  from  experimental  data  and  it  is 
implicitly  assumed  that 
J 


o 

u(-  ) 

0  ^  t^  t . 


^deploy)  +  j 


(swing) 


r(stance) 


0<tl<t2<tf  . 


u(-  ) 
0^  t<  t 


u(') 


t  ^  t  5;  t 
C1  ‘2 


u() 

t^t^tf 


In  general,  the  optimal  control  strategy  over  the  entire  interval 
[0,t  ]  is  not  equal  to  the  concatenation  of  the  optimal  strategies  over 
the  successive  sub-intervals.  However,  it  appears  that  the  very 
special  nature  of  the  locomotion  problem  renders  the  sub-optimal 
solution  close  to  the  optimal  one.  In  both  the  deploy  and  stance  por¬ 
tions,  the  contributions  due  to  the  ankle  moment  and  reactions  are 
the  major  determinants  in  the  control  moment  u(-  )  calculation.  In 
addition,  the  moments  in  the  stance  portion  are  determined  by  an 
algebraic  method  without  an  optimization  procedure  and  so  the 
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contribution  j(stanCt  )  is  essentially  constant,  irrespective  of 

optimization.  The  cost  j(dcploy)  a  sma]q  fraction  of  the  total 

1 1 

cost  because  the  subinterval  [0, t^]  is  short  (i.e.  - —  ^  0  1  to  0.  1  5) 


and  the  control  moments  are  also  influenced  by  the  ankle  moment  and 


(swing) 


reactions.  Only  the  term  J 


can  be  significantly  minimized 


via  optimization.  Thus,  it  appears  that  the  sub-optimal  solution 
should  be  very  close  to  the  true  optimal. 

3.  Necessary  Conditions  of  Optimality;  Method  of  Solution. 

With  the  problem  formulation  completed,  the  next  step  is  to 


invoke  the  necessary  conditions  of  optimality.  Solution  for  the  deploy 


and  swing  portions  are  first  considered  as  they  are  of  a  similar 
nature.  The  technique  is  to  convert  the  original  constrained  optimi¬ 
zation  problem  into  a  sequential  unconstrained  problem.  Let  S(x;t) 
denote  the  state  constraint  (equality  or  inequality)  in  either  portion. 
Define  an  extra  state  variable  x^(t)  such  that  it  satisfies  the  differential 
equation 


*5(t)  =  sgn(S)  •  SZ(x;  t) 


(10) 


1  if  S  =  0 


sgn(S)  =  ' 


(11) 


0  if  S  ^  0 


for  O^t^t^;  t^^t^t^.  The  original  state  vector  satisfies,  of 
course,  the  dynamic  equation 


(12) 


With  the  initial  condition  set  to  zero,  i.e.  x^(t^)  =  0  for  swing 
x^(0)  =  0  for  deploy  respectively,  then 
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(13) 


The  state  constraint  is  thus  satisfied  A  mathematical  proof  of  the 
above  penalty  function  technique  is  given  by  Lele  and  Jacobson  [44]. 
As  far  as  the  initial  and  terminal  conditions  are  concerned,  a 
quadratic  penalty  function  is  used.  With  these,  the  necessary  condi¬ 
tions  of  optimality  can  be  separately  derived  for  the  deploy  and  swing 
portions . 

In  deploy,  we  have  the  modified  performance  criterion  for  the 
sequential  unconstrained  problem,  i.  e. 


A 


Min  J  ,  =  cr 
d 


u 


t 


+  (x4(tx)  -  d4)2}  +  ~ 


(14) 


0 


Specifically,  we  solve  (14)  for  a  monotonically  increasing  sequence 
0<  cr  ,  <  cr  ?  < - <  cr  an<3  subject  to  the  equations 


(15) 


d 


The  signum  function  Sgn(S  )  =  1  because  the  state  equality  constraint 
is  effective  throughout  the  interval.  Define  a  variational  Hamiltonian 


(16) 


The  X's  are  the  influence  functions  or  the  adjoint  variables  of  the 
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optimization  problem.  The  extremal  condition  is 

4 


9H 

9u. 


8f. 

=  0  =^r.u  +  )  X  — — =  0  for  i  1,2 

ii  ^  j  3u 

j  =  l 


(17) 


rlUl  +  6  (A2X3  +  SV  =  ° 


r2U2  +  5  (t7X3  +  t8X4)  '  ° 


(17a) 


Note: 


2  d 

a  H 


8u. 


tc—  -  r.  >0  for  r.  >  0:  i  =  1.2.  The  influence  equations 

2  i  i  7  ? 


and  their  terminal  conditions  at  t^  are 

r  -x.  -  aH<1 


i  9x. 


V  af;  a  p>cd 

=  2  V-  +2V  S  ‘ 


=  1,— ,4 


i=l 


-X  =  0 

D 


9^ 


dfi .  9^ 

X2<*l>  =  2<rk<*l  —  +* 


2  2  8x2 


^ 3 (l  1 )  =  2ak(x3(tl)  '  d3) 


^4^T  _  2<r  )  ~  da^ 


k'  4'  1‘ 


v  S(ti>  "k 


9S 


Explicitly,  expressions  for  the  X-equations  and 


\ 


(18) 


are  given  as  follows 
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X1  -5  {(A2R3i +t7R4i)^3  +  <t7R3i  +t8R4i)K4)  +  2X5S  '  9x 


^2  TT  1  A2R32  "C3l5^R4  tu2^  +t7  ’  R42  ~5~  ■  A2^R  3 +U1  t7^R4+ U2^ 

X4 

+  -5-  ’  (t7R32  '  C3t5(R3  +ui}  +t8R42  '  ZC3t5(R4+  u2} 

+  •  [t7(R3+Ul)  +tg  •  (R4  +  u2)]}  +  2^5Sd-  |~ 

-  X3  =  \l  +  -5  {(A2X3  +  t_,X4)  •  R33  +(1-7X3  +t8^4)  '  R4  3  } 

R 

"  X4  =  x2  +  - {  a2x^  +t?x4) 

-  X5  =  0 

•  « 

R  3  =  _C3t5x4(x4  '  2x3>  "  (C111  '  C2t3^8  +  g) 


R  4  =  -  C3t5x3  _  C3t6(g  +  g) 


R 3 1  "  (Clt2  +  C2t4^t10  ;t10  =  g+^ 


R41  CZt4t10  +  VZY  "  VlX 


R3Z  =  "  C3t6x4(x4  '  Zx3}  +  C2t4  '  *10 


R 


4Z 


C3t6x3  CZt4t10 


R33  =  2C3t5x4 


R43  _2C3t5X3 


R 34  ~  _2C3t5^x4  "  x3^ 


Z  Z  Z 

S  1  +S2  '  d 
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S1  =  el  -  ljt2  -  l2t4  -  g(t) 


S-,  =  e,  -  £,tn  +  -  v  t 

2  2  11  Z  3  o 


8xx  "  2Si  fci  '  £2t3)  ’  2S2^1t2  +  £2t4) 


|f-=2i2(Slt3  +  S2t4) 


6  "  A2t8  "  l7 


6  #  |§-  =  2C2t,t 

2  3  5  6 

Similarly,  in  the  swing;  portion,  we  have  the  augmented  performance 
crite  rion 


4  2  f2 

MinJs  =  o'  k(x5(t2)  +_  (Xj(t2)  -  d.)2}  +  vA 

u  i=l  J 


dt 


x2(t) 


1  \  ,  2  2. 

+  2  )  (rlul  +  r2u2)dt 

\ 


(19) 


forO  <  c r  ^  < - cr  — ►  co  and  >  Y;?  - >  y^— ^  0_^  .  The  dynamic 

equations  are  now 

x(t)  =  f(x,  u;  t)  ;  x(t^)  s.t.  ^  ^  =  0,  ^  =  ^  from  deploy 


x5(t)  =  Sgn(SS)  •  (Ss(x;t))  ;  x^tj)  =  0. 


(20) 


The  presence  of  the  "soft1’  penalty  is  to  ensure  that  x^(t)  >  0  during 
the  entire  swing  portion.  This  requires  that  the  initial  and  all  iterative 
trajectories  be  feasible.  Similar  to  the  deploy  portion,  define  a 


variational  Hamiltonian  as 
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H  (x,u,X;cr  yk;t)  +  rZuZ^  +  x  (t)  +  l^Vj1j^X,u;t) 


j=» 

+  X5(t)  •  Sgn(SS)  •  (Ss(x;t))<i 


(21) 


Again,  the  extremal  condition  implies 

4 

rS  i 


9f. 


(H  )u-0=^>r1u1  ^  -0;  i  =  1,2 


(22) 


j  =  l 

The  equation  set  (17a)  is  thus  obtained.  The  influence  equations  and 
their  appropriate  boundary  conditions  are 


.  V  a£i  SB  9Ss  YV -6(2,0  ^ 

-\=ZXj  KL+2X5Sgn(S  »  '  S 
j=l 


3x.  2 

1  x^(t) 


-x5  =  0 


\(t2)  =  2tr  k(xi(t2)  '  di) 


(23) 


X5(t2)  "  °  k  J 

for  i  =  1, - ,4  and  6  (2,  i)  -  1  if  i  =  2  and  0  otherwise.  The  expressions 

for  the  X-equations  are  similar  to  those  for  the  deploy  portion  except  for 
the  substitution  of  inequality  for  equality  constraint  expressions. 
Expressions  related  to  the  constraint  are 

-ijtj  +£zt3  -  d  sin  (3(t) 


< 


9S 

9x^ 


=  +  d  sin  P(t) 


(24) 


|§!,0;  flsO 

9x^  9x^ 
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For  both  the  deploy  and  swing  portions,  we  thus  have  a  two-point 

boundary  value  problem  (TPBVP)  -  (15,  17a,  18)  and  (20,22,23) 

for  the  respective  unconstrained  optimization  problems. 

Solution  of  the  stance  portion  differs  from  the  other  two  stages. 
The  leg  in  stance  carries  the  body  weight  while  the  inertia  of  the 
trunk  and  forward  motion  of  the  swinging  leg  carry  the  entire  body 
forward.  The  result  is  a  ’'sinusoidal 11  trajectory  for  the  hip  and  the 
foot  is  almost  stationary  during  the  roll-over.  With  prescribed 
motion  for  the  hip  and  ankle,  motion  of  the  thigh  and  shank  are 
determined  through  these  kinematic  constraints.  The  result  is  that 
the  effective  hip  and  knee  moments  can  be  found  without  the  need  of 
minimization  of  the  performance  criterion;  an  algebraic  method  is 
sufficient  for  the  solution  in  this  portion.  The  equality  constraints 
on  the  ankle  motion  are 


sfS(x;t)  =  g(t)  +  +  £~t  +  p(t)  -  e  =  0 


1 


V.  S^b(x:t)  =  v  (t  +  t  )  +  tx  t,  -  1  -  q(t)  =  0 

2  7  o  o  11  2  3 

for  t^  ^  t  ^  t^.  Differentiating  once  with  respect  to  t,  we  have 
Sp  (X- 1)  -  gU)  +  p(t)  +  z11x3  +  z12x4(t)  -  0 

L  ^2  S  (x»  t)  -  v0  '  4(t)  +  Z21X3^^  +  z22x4^  “  0 


12  2  4 


(5  ) 


(25) 


where 


z  = 


Z11  z  1  2 


Z21  Z22 


l-l  t  4-  1  t 

11  2  3 

"  lZll \ 

at  +  a  t 
l12  T  2  4 

-  vJ 

The  equation  sets  (5')  and  (25)  form  a  decoupled  system  of  four 
equations  in  the  state  variables  (x^, - ,  x^)  at  each  instant  of  time  in 


64 


i 

the  stance  interval.  Note  that  (5  )  is  a  nonlinear  equation  set  in  x^ 

and  x~»  while  (2  5)  is  linear  in  and  x  .  with  the  coefficients  z 
^  3  4  ij 

(i,  j  -  1,2)  as  functions  of  x^  and  Differentiating  once  more  with 

respect  to  time, 

h 


/  *  *  r  s  \ 

si  ) 

sV  1 

z  / 

g  +  p 


-  q 


+  z 


V2 


L 


\ 


=  0 


y(x4-x3)  j 


(26) 


where 


£2t4' 


’ £2t3/ 


Now,  x^  and  x^  contain  the  actuating  moments  u^  and  u,,  explicitly 
from  the  original  dynamic  equations 


or 


X1  "  x3 


x2  “  x4 


=  —  {A^R^  +  u^  +  M^  +  ^\l\  ”  £2t3^Y  +^lt2  +  £2t4^  +t7^R4 


+  u2-Ma  +  i2t3Y-i2t4X)} 


(27) 


x4  =  '{t^Rj+Uj  +Ma  +  (^i1!  "  *2*3^  +^’1t2  +  ^2t4^  +  t8^R4  +  u2 


-  M  +  £0t,,Y 
a  2  3 

’  £2 

}K\ 

L 

+  A  " 1 

\  R4  i 

\U2  J 

iA. 


A'1  -  I 
’  A  "6 


yi  h 


=  A  t,  -  t_ 
2  8  7 


(271) 


Combining  (26)  and  (27'),  the  moments  can  be  expressed  as: 
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u  =  -  A 


(  2 


x. 


..  V  A 

-  (g  +  P)\ 


+  AZ 


-1 


L 


Jx4  "  x3) 


1 


q(t) 


V 


/ 


(28) 


and 


lk  3  \  ( 

S' 

x3  \ 

I  * 

1  -  (g  +  p) 

hz'1 

L 

!  ^ 

( 

\ij 

V(X4'X3’ ) 

1 

V  5(t)/ 

(29) 


To  generate  the  trajectories  in  the  stance  portion,  we  proceed 


forward  from  t  =  t^  to  t  =  t^  in  small  positive  increments,  A.  At  each 


instant,  solve  equation  set  (51)  for  x^  and  numerically  (this  is 
essentially  finding  the  roots  of  nonlinear  equations).  With  the  x^  and 
x2  values  determined,  the  matrix  Z  is  known  and  so  the  variables  x^ 


and  x.  can  be  computed  from  (25).  Equation  (28)  then  generates  the 


hip  and  knee  moment  profiles. 
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E.  Numerical  Simulation  and  Results 


1 .  Parameters  and  Pertinent  Data 

Having  formulated  an  optimal  programming  theory  for  biped 
locomotion,  the  next  step  is  to  perform  a  numerical  simulation  of  the 
gait  based  on  the  proposed  theory.  This  is  a  crucial  step  because 
only  the  results  can  demonstrate  the  validity  and  relevance  of  the 
theory.  In  general,  the  walking  pattern  of  each  individual  is  a  little 
different  from  that  of  another  person  because  of  the  skeletal  structure 
and  physiological  conditions.  However,  there  are  characteristic, 
qualitative  features  which  are  observable  in  the  human  gait  for  non- 
pathological  cases.  The  objective  of  the  numerical  experiment  is  to 
simulate  these  common  features, 
a .  Physical  Parameters. 

The  values  of  the  parameters  used  in  the  present  computation 
are  based  on  those  used  by  Galiana  and  Milsum  [33]  in  their  EMG- 
experiment.  Typical  values  for  the  mass,  length  and  moment  of 
inertia  for  the  two  segments  of  the  lower  extremity  are  listed  as 
shown  below.  It  is  noted  that  their  values  are  derived  from  the 
commonly  used  !,pe rcentage 11  values  of  Fisher  and  Dempster  in 
Lissner  [32].  Values  of  the  composite  parameters  used  in  the  mathe¬ 
matical  model  are  also  listed 

Parameter  Value  Remarks 


Body  Weight,  W  130  pounds 

4.  05  slugs 


Mass,  A 


Aq  =  W/g 


Mass  of  Thigh,  M^  0.  389  slugs 
Mass  of  Shank,  M^  0.  2595  slugs 


9 .  6%  of  A 


6.  4%  of  A  (  includes  mass  of 
°  foot) 
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Parameter 

Value  Remarks 

Moment  of  Inertia  of  Thigh,  1^ 

1 . 4  x  10  ^slugs 

Moment  of  Inertia  of  Shank,  I^ 

-2 

1.882x10  slugs  (shank-foot  combination) 

Link  Length  of  Thigh, 

1.  3524282  ft. 

Link  Length  of  Shank, 

1.  2618726  ft. 

Distance  from  Hip  Joint 

0.  586  ft.  Dempster’s  average  data 

to  C.  G.  of  Thigh,  a 


Distance  from  Knee  Joint 

0.  547  ft.  Dempster's  data 

to  C.  G.  of  Shank,  a  ^ 

Length  of  Foot 

/*) 

0.  55  ft.  Effective  length  in  deploy 

0.  65  ft.  Effective  length  in  swing 

Composite  Parameters  in  Simulation 

=  1^  +  m^a  ^  5993547648  slug-ft. 


.  t  i  ^ 

A2  =  I2  +  m 2a  2 

=  0.  1414525055  slug-ft.  2 

Ci  =  m i  4  ™ 2} \ 

=  0.  564730024  slug-ft.  2 

C  2  ~  mZa  2 

=  0.  1735743239  slug-ft.  2 

C3  =  C2  £1 

=  0.  2347468104  slug-ft.  2 

b.  Initial  and  Terminal  Conditions. 

The  starting  values  for  the  angular  displacements  and  velocities 
(i.  e.  initial  states)  are  derived  from  Beckett’s  [28]  numerical  experi¬ 
ment  These  are 


(*)  During  deploy,  the  foot  rotates  about  its  ball.  The  distance  between 
the  ankle  joint  and  center  of  rotation  is  shorter  than  the  foot-length 
during  swing  when  the  toes  are  also  taken  into  consideration. 
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x  (0)  =  0.  1  1  170107  radians  (6.4°) 

x^(0)  -  0.  22689^8  radians  (1  3°) 

x^(0)  =  x^(0)  =  0.  radians/sec. 

x^(0)  =  x^(0)  =  1  .  68  radians/sec. 

At  the  end  of  the  double-step,  the  terminal  values  are  chosen  as 
x^(2.  0)  =  0.  799  radians  (46°) 

x2(2.  0)  -  0.  052  radians  (3°) 

x3^’  ^  “  0.0  radians/sec. 

x^(2.  0)  -  0.  0  radians/sec. 

The  initial,  intermediate  and  terminal  times  for  the  three  portions 

are  set  at  t^  -  0;  t^  0.  2  sec  (deploy);  =  1.0  sec  (swing);  =  2.  0 

(stance).  Here,  we  assume  a  10%  interval  for  deploy  and  restraint 

activity.  The  swing  (or  support)  occupies  40%.  Variation  in  the 

temporal  patterns  modify  the  numerical  results  but  not  the  methodology 

c .  Prescribed  Trajectories. 

For  the  hip  motion,  we  have 

f(t)  =  v  (t  4  t  )  feet 
o  o 

v  -  2.275  ft .  /sec . 
o 

t  =  0.  24698953  sec. 
o 

g"(t)  =  e^  ■  s^n  (t  +  P  t^)  feet 

e  =  2.  73581  72  feet 
o 

t^  =  2.  0  sec. 


P  =  0.  525 
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e,  =  e  =  2.  73  58172  feet 

1  o 

e,  =  -  v  t  =  -  0.  399078  feet 

2  o  o 

To  describe  the  almost  stationary  behavior  of  the  foot  in  stance 
portion,  we  rely  on  the  experimental  profiles  by  Bresler  and  Fran^el 
[13],  and  other  experimenters  [10,  12],  (Fig.  11).  For  the  present 
simulation,  simple  curve  fitting  on  the  experimental  data  is  used  as 
our  important  objective  is  to  illustrate  our  approach  of  study.  Better 
results  could  be  obtained  through  interpolation  or  spline  function 
approximation.  In  the  stance  portion  simulation,  the  following  prescribed 
trajectories  are  used  for  the  ankle  motion: 

q(t)  =  2.00+  0.35(l-e'6-tlT)  -  T=t-1.0;1.0^t«2.  0. 

0.17  +  0.15Sgn(r7)*r7^  1.0^t^l.3 


P(t) 


0.  17 

0.  17  +  0.  25(t-  1 .  0)  +49.  1  •  [0.  25(t-l)4  -  0.  755(t-l)3 


1.  3<t^  1.  55 


+  0.  85(t-l)  -  0.  42 1  (t-  1)] 

(/?  =  t  -  0.  3  ;  Sgn (r?)  =  { q 


1.  55«t  -'2.0 


n  §  o 

r?  >  0 


The  particular  choice  of  (p,  q) -profiles  is  based  on  the  form  of  the 
experimental  curves  as  well  as  continuity  at  the  matching  points, 
d.  Pressure  Transfer  Curve. 

Equation  (46)  of  mathematical  modelling  proposes  a  linear 
expression  for  the  transfer  characteristics  of  the  pressure  center  inthe 
stance  portion.  The  distance  travelled  from  heel  strike  to  rotation  about 
the  ball  of  the  foot  is  equal  to  the  effective  foot  length  in  deploy,  i.  e. 
d  =  0.  55  \  In  approximately  mid-stance,  the  ankle  joint  is  directly 
above  the  pressure  center.  These  two  pieces  of  information  give  us 
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the  required  expression: 

x  (t)  =  2.  35  +  0.  9  (t  -  1 . 40)  feet  1 . 0  ^  t  ^  2.  0 

2.  Numerical  Computations. 

From  the  conditions  of  optimality,  a  first  order  algorithm  is 
employed  for  iterative,  numerical  solution.  The  algorithm  is  coded  into 
Fortran  IV  programming  language  to  run  on  the  IBM  360/65  digital 
computer.  The  program  possesses  the  following  features: 

a.  The  program  has  two  options: 

(1)  a  one -dimensional  minimization  subroutine  (QUADCL)  used  in 
conjunction  with  conjugate  gradient  or  steepest  descent  algorithms. 

(2)  a  simple  gradient  procedure  with  a  pre  -  specified  improvement 
step- size . 

b.  A  fourth  order  Runge-Kutta  method  with  fixed  integration  step- 
size  is  used  for  integrating  the  state  and  adjoint  equations.  In  each 
integration  step,  the  control  histories  u(t),  state  trajectories  x(t)  and 
ground  reactions  (X(t),  Y (t ) )  are  represented  by  their  values  at  the  end 
points.  Presumably,  such  a  ’’discretization"  could  have  effect  on  the 
quality  and  convergence  rate  of  the  computation. 

c.  To  satisfy  the  state  constraints,  the  computation  is  done  in 
series  with  increasing  positive  values  of  the  weighing  factor,  <r  . 

d.  The  computation  is  carried  out  in  singie  precision  arithmetic. 

The  numerical  computation  is  involved  because  the  dynamic 

equations  involve  five  state  variables  and  they  are  highly  non-linear. 
More  efficient  algorithms  such  as  second  order  optimization  methods 
are  very  complicated  to  program.  In  the  literature,  little  has  been 
reported  on  the  computational  experience  involving  problems  of  such 
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complexity.  Besides,  the  primary  purpose  of  the  computation  is  to 
establish  the  relevance  of  the  theory  to  biped  locomotion  and  simulate 
the  common  features  of  gait.  These  considerations  prompt  the  use 
of  the  first  order  algorithm  for  an  adequate,  approximate  solution  to 
the  problem.  For  each  value  of  <r  ,  the  computation  involves  a  large 
number  of  iterations  before  a  satisfactory  solution  is  obtained.  The 
procedure  with  the  QUADC L.- search  option  involves  more  computing 
time  per  iteration  in  general  than  the  simple  fixed-step  gradient 
procedure.  However,  these  two  options  were  used  more  or  less 
alternatively  to  achieve  improvement  (cost  reduction)  because  of  the 

non-linearity  of  the  problem.  An  acceptable  solution  is  obtained  with 

3  2 

<r  =10  for  the  deploy  portion  and  cr  =  10  for  the  swing.  An 

interesting  observation  is  that  the  swing  portion  with  an  inequality 

constraint  is  comparatively  simpler  to  compute  than  an  equality 

constraint  as  in  the  deploy  portion. 

3 .  Results  and  Discussion. 

Results  of  the  simulation  are  summerized  in  the  various  graphs. 
To  facilitate  discussion,  comparison  is  made  with  the  often-quoted 
experimental  works  of  the  Ebe rhart-Inman  group  at  Berkeley  [10,  12, 
13].  In  general,  the  results  can  be  classified  into  three  groups  for 
discus  sion. 

a.  Angular  Displacements  and  Related  Results. 

Results  pertaining  to  this  part  are  depicted  in  the  graphs  No.  1 
to  No.  7.  The  state  trajectories  x^(t)  and  x>(t)  together  with  their 
time  derivatives  portray  the  angular  motion  of  the  thigh  and  shank  in 
the  sagittal  plane.  Note  that  the  curves  are  so  plotted  that  they  begin 
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at  heel  strike  (i.  e.  stance  portion)  for  comparison  purposes. 

The  x^(t) -trajectory  describes  the  motion  of  the  shank  with 
reference  to  the  thigh.  At  heel  strike,  the  knee  joint  is  extended 
so  that  the  x7(t)  angle  is  at  a  small  positive  value.  During  the  brief 
restraining  portion  of  stance,  the  angle  x^  increases  so  that  this  flexing 
of  the  knee  absorbs  the  impact  of  body  weight  passing  on  to  the  leg  in 
stance.  The  angle  then  decreases  until  deploy  when  it  starts  to  increase 
rapidly.  The  behavior  that  the  knee  is  first  flexed  at  heel  strike  and 
then  flexed  again  during  deploy  is  known  as  f,knee  locking”  -  a  fact  well 
documented  in  experimental  studies.  In  the  swing  portion,  x^(t)  reaches 
a  peak  value  of  roughly  80°  and  then  decreases  rapidly  to  almost  zero 
in  0.  5  second.  The  behavior  resembles  the  ballistic  action  of  a 
pendulum  as  illustrated  clearly  by  the  velocity  profile  x^(t)  in  graph 
No.  2.  In  the  pioneering  works  of  Fischer,  it  was  conjectured  that  the 
action  of  the  leg  in  swing  is  purely  due  to  gravity,  but  this  was  disputed 
by  others  later  in  connection  with  muscle  activity  studies.  From  the 
present  optimization  standpoint,  it  appears  that  the  effective  moments 
cooperate  with  gravity  for  lesser  expenditure  of  work  while  still 
satisfying  the  state  constraints.  From  the  graph  No.  17,  it  is  apparent 
that  the  constraint  is  well  satisfied.  The  state  trajectories  (i.e. 
angular  displacements  and  velocities)  are  well  matched  for  the  multi¬ 
arc  programming  problem.  This  is  achieved  by  using  the  terminal 
conditions  of  one  stage  as  the  initial  values  for  the  following  in  actual 
computation.  However,  there  exist  small  discontinuities  for  the 
acceleration  and  moment  profiles  at  the  stance/deploy  and  deploy/swing 
junctions.  The  reason  for  the  mismatch  is  primarily  numerical  as  the 
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penalty  technique  used  in  the  computation  seldom  satisfies  the  desired 
conditions  exactly.  For  the  purpose  of  illustrating  the  theory,  the 
"numerical"  discrepancy  is  not  critical.  Similarly,  another  interesting 
point  is  that  the  ~accele ration  curve  (graph  No.  3)  exhibits  a  "peaky" 

behavior  in  mid-stance.  The  probable  reason  is  that  in  our  algebraic 
method  of  solution,  a  piecewise  polynomial  representation  is  used  for 
the  ankle  trajectories.  Another  reason  for  the  slight  mis-match  of 
acceleration  profiles  could  be  that  the  overall  optimization  problem  is 
decomposed  into  a  number  of  optimization  problems  over  sub-intervals 
defined  by  the  different  phasic  periods.  This  results  in  an  overall 
sub-optimal  solution  unless  the  joins  of  the  sub-intervals  happen  to 
occur  at  exactly  optimal  locations  (which  is  unlikely  since  these  are 
experimentally  obtained  data)  and  particular  care  is  taken  in  the 
matching  of  the  adjoint  variables. 

The  hip  trajectory  x-^(t)  portrays  the  familiar  pattern  as  observed 
in  experimental  studies.  In  the  stance  portion,  the  angle  x^(t)  decreases 
more  or  less  uniformly  while  the  hip  joint  is  describing  a  sinusoidal 
path.  This  characterizes  the  roll-over  action  of  the  upper  part  of  the 
body  (HAT- section)  about  the  weight-bearing  leg.  During  deploy  and 
subsequent  swing,  the  angle  x^(t)  shoots  up  to  a  maximum  value  and  then 
coasts  until  heel  strike.  The  combined  motion  of  the  thigh  and  shank 
clears  the  toe  sufficiently  above  the  ground  for  swing  motion.  As  for 
the  x^(t) -trajectories,  the  x^(t)-trajectories  are  matched  at  the  junction 
points  except  for  the  acceleration  curve.  The  reason  is  again  numerical 
as  discussed  above. 

The  angle  a  (t)  (Fig.  12)  describes  the  rotation  of  the  foot  about 
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its  ball  during  deploy  action.  Graph  No.  4  shows  the  result  of 
computation  based  on  the  equality  constraint  relationship.  The  angle 
a  (t)  has  a  monotonically  increasing  profile,  reaching  a  limiting  value 
of  about  1.  2  radians  (~  68°)  at  the  end  of  deploy.  A  value  of  70°  was 
used  for  specifying  the  terminal  deploy  conditions.  The  velocity  and 
acceleration  curves  likewise  show  a  monotonically  increasing  trend, 
but  the  acceleration  trajectory  appears  to  be  more  irregular.  This  is 
because  the  state  constraint  is  satisfied  approximately  as  shown  in 
the  residue  error  curve  in  graph  No.  16. 

The  angular  displacements,  velocities  and  accelerations  con¬ 
stitute  the  main  results  of  the  computation  upon  which  related  quantities 
can  be  calculated.  Graph  No.  5  shows  a  phase  plane  diagram  for  the 
angular  displacements  x^(t)  versus  x^(t)  during  a  double-step.  In  the 
paper  by  Marie  and  Vukobratovic,  a  special  type  of  compass  gait  was 
considered  in  which  the  hip  motion  is  stationary  in  the  vertical  direction 
and  the  foot  is  in  continuous  contact  with  the  ground.  The  resulting 
realization  of  the  gait  pattern  is  that  the  thigh  and  shank  angles  are 
sinusoidal  in  a  double-step.  This  gives  an  elliptic  phase  plane  plot 
with  the  points  equally  spaced.  The  present  results  simulate  much  more 
realistically  the  biped  gait  and  the  plot  indicates  deformation  of  the 
curve  in  the  stance  portion.  Note  that  the  points  are  not  uniformly 
spaced,  showing  different  rates  of  activity  in  the  three  portions. 

The  horizontal  and  vertical  displacements  of  the  joints  are 
plotted  in  graphs  No.  6  and  No.  7.  A  double-step  takes  two  seconds 
during  which  there  is  an  overall  horizontal  displacement  of  4.  55  feet. 
This  gives  a  walking  speed  of  about  1.  6  miles  per  hour,  a  case  of  slow 
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level  walking.  The  hip  moves  at  a  uniform  rate  while  the  knee  and 
ankle  move  through  4  feet  during  the  swing  portion.  Graph  No.  7 
illustrates  the  vertical  displacement  of  the  joints.  A  double  sinusoid 
with  0.  9  inches  in  amplitude  is  used  to  simulate  the  hip  motion.  In  a 
more  accurate  description,  a  Fourier  series  representation  can  be  used 
to  fit  the  experimental  profile.  In  fact,  the  amplitude  and  phase  of  the 
various  components  have  been  studied  to  reveal  the  abnormal  properties 
of  the  pathological  gait.  In  the  normal  case,  the  double  sinusoid  is  the 
principal  component.  During  stance,  the  knee  joint  is  almost  stationary, 
so  is  the  ankle  trajectory  as  prescribed  by  the  p(t) -profile .  The 
displacement  at  initial  stance  is  prevailingly  restraining  action  to 
absorb  impact  and  the  gradual  increase  towards  deploy  implies  the 
transfer  of  the  pressure  center.  In  subsequent  deploy  and  swing  motion, 
the  vertical  motion  of  the  knee  and  ankle  is  comparatively  larger.  Note 
that  the  overall  vertical  displacement  is  almost  zero  in  a  double  -  step, 
b.  Reaction  Forces  and  Moments  at  the  Joints. 

The  foregoing  section  describes  the  kinematic  behavior  of  the 
lower  extremity  in  level  walking.  To  complete  the  description, 
dynamic  quantities  such  as  the  forces  and  moments  at  the  various  joints 
must  be  determined.  It  is  only  with  the  dynamic  information  that  a 
realistic  description  of  biped  locomotion  can  be  said  to  have  been 
achieved.  In  the  weight  bearing  portions,  the  foot  is  subjected  to  very 
large  ground  reaction  forces  due  to  the  weight  and  acceleration  of  the 
upper  part  of  the  body  (HAT-section).  With  these  large  forces 
appearing  at  the  inputs  to  the  dynamic  equations,  it  is  not  difficult  to 
realize  that  they  are  the  dominant  quantities  in  determining  the 
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actuating  moments  and  reaction  forces  at  the  various  joints.  In  passing, 
we  note  that  the  numerical  integration  of  differential  equations  having 
inputs  of  a  large  rapidly  varying  magnitude  is  not  trivial. 

The  present  model  assumes  the  ground  reactions  as  known  inputs 
to  the  dynamic  equations.  This  condition  amounts  to  prescribing  the 
behavior  of  the  HAT-section  and  then  investigating  the  walking  patterns 
of  the  lower  extremities  using  optimization  theory.  It  is  certainly 
consistent  with  the  concept  of  decomposing  the  human  postural  study  into 
the  stability  problem  of  the  HAT-section  and  the  locomotion  activity 
study  of  the  two  legs,  i.  e.  control  through  a  hierarchy  structure.  In  an 
analysis  using  a  complete  high  order  model,  the  ground  reactions  are 
determined  from  the  resulting  motion  of  the  various  segments.  For  the 
present  study,  the  experimental  profiles  by  Inman  [31]  (normalized  to 
body  weight)  are  used  as  shown  in  graphs  No.  8  and  No.  9.  Note  that 
the  vertical  reaction  has  a  double-peak  behavior  in  the  stance -deploy 
portions.  The  first  occurs  after  heel  strike  when  the  body  rolls  over 
the  weight-bearing  leg.  The  second  peak  occurs  in  the  deploy  portion 
when  the  leg  pushes  the  body  up  in  preparation  for  the  swing  motion 
In  both  graphs,  it  is  apparent  that  the  reactions  at  the  ankle,  knee  and 
hip  joints  are  close  to  the  ground  reactions.  The  differences  are  due 
to  the  inertia  and  mass  of  the  thigh  and  shank  segments.  The  ground 
forces  are  zero  in  the  swing  portion. 

Results  pertaining  to  the  moment  quantities  are  presented  in 
the  graphs  No.  1  0  to  No.  12.  The  ankle  moment  has  the  greatest 
values  in  the  weight  bearing  portions.  The  manner  in  which  the  ankle 
moment  is  calculated  is  approximate.  Based  on  a  "quasi- static  " 
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assumption,  the  inertial  and  gravitational  effects  of  the  foot  are 
neglected  as  they  are  1%  or  less  of  the  contribution  of  the  ground 
reactions.  In  addition,  a  linear  function  is  used  to  describe  the 
transfer  of  pressure  center  in  stance.  Despite  the  crudeness  of  the 
calculation,  the  ankle  moment  agrees  qualitatively  well  with  the 
experimental  results  of  Bresler  and  Frankel  [13].  This  is  encouraging 
because  the  approximate  calculation  avoids  the  complexity  of  the  high 
order  model.  With  better  representation  (e.  g.  spline  polynomial 
interpolation)  of  the  prescribed  trajectories,  better  results  are  antici¬ 
pated.  The  ankle  moment  and  ground  reaction  contributions  are  the 
dominant  quantities  in  determining  the  hip  and  knee  moments  in  the 
stance  and  deploy  portions  (graphs  No.  11  and  No.  12).  This  is  very 
interesting  because  the  inertial  and  gravitational  terms  are  the  most 
complicated  to  calculate  and  yet  they  are  only  second  in  importance  in 
the  moment  calculation.  Thus,  the  use  of  the  simplified  decoupled 
model  together  with  kinematic  constraints  is  not  as  crude  as  it  may 
initially  appear.  However,  neglecting  the  ankle  moment  and  ground 
reactions  makes  a  big  difference  in  the  moment  calculation.  This, 
then,  should  explain  the  shortcomings  of  the  papers  by  Beckett, 

Wallach,  Galiana  and  Vukobratovic  et  al. 

The  hip  and  knee  moments  are  also  plotted  in  graph  No.  10. 

The  greatest  activity  occurs  during  stance  following  heel  strike  and 
at  deploy.  Activity  in  the  swing  portion  is  addressed  to  maintaining 
the  leg  above  the  ground  (i.e.  to  satisfy  the  state  inequality  constraints) 
while  utilizing  gravity  as  much  as  possible). 

The  knee  moment  profile  shows  a  typical  pattern  associated  with 
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the  "knee  "locking "  behavior  discussed  earlier  in  a.  The  sign 
convention  for  the  moment  is  such  that,  when  negative,  the  forces 
pass  behind  the  knee,  tending  to  flex  it;  while  a  positive  moment 
tends  to  hold  the  knee  in  a  "locked"  position.  The  present  profile 
agrees  fairly  well  with  experiment  [13].  At  heel  strike,  the  leg  is 
quickly  stabilized,  then  unlocked,  and  once  again  locked  for  the 
deploy  portion.  Bresler  [13]  pointed  out  that  the  knee-locking  behavior 
allows  one  to  move  forward  with  a  minimum  raise  of  center  of  gravity, 
and  therefore  less  expenditure  of  work. 

Like  the  ankle  and  knee  moments,  the  hip  moment  again  agrees 
with  the  experimental  pattern.  The  negative  moment  following  heel 
strike  is  responsible  for  the  roll-over  of  the  body,  indicated  by  a  mono¬ 
tonic  decrease  of  the  x^(t)-ankle.  Toward  the  deploy,  the  hip  angle 
starts  to  increase.  The  greatest  values  at  initial  stance  and  deploy 
are  opposite  in  sign.  Based  on  the  idea  of  two  basic  configurations  for 
the  biped  gait,  the  two  portions  of  "maximum"  activity  correspond  to  the 
"restraint/deploy"  configuration  (Fig.  6).  From  the  dynamic  equation 
for  the  HAT-section,  the  trunk  of  the  body  is  subjected  to  the  action  of 
the  hip  moments.  But  the  interesting  fact  that  the  hip  moments  of  the 
two  legs  are  more  or  less  of  the  same  order  of  magnitude  but  opposite 
in  sign  leads  to  the  concept  of  moment  annihilation  (i.  e.  cancelling)  of 
the  hip  moment’s  effect  on  the  HAT-section.  This  consequently  requires 
far  less  stabilizing  moment  effort  for  maintaining  up- right  posture.1  In 
the  same  vein,  the  "complementary v  profiles  for  the  hip  reaction 
forces  (graphs  No.  8  and  No.  9)  create  a  couple  action  in  the  frontal 


j  i.e.  ,  the  profiles  at  deploy  are  reflections  of  those  at  initial  stance. 
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plane,  facilitating  the  transfer  of  body  weight  from  the  deploying  leg 
to  the  restraining  one  following  heel  strike, 
c.  "Quality11  of  Numerical  Computation. 

As  stated  earlier,  the  important  objective  of  the  numerical 
experiment  is  to  demonstrate  the  validity  and  relevance  of  the  optimal 
programming  approach  to  the  study  of  biped  locomotion  The  require¬ 
ment  o £  an  adequate  solution  together  with  complexity  of  programming 
prompted  us  to  use  a  first  order  algorithm  and  appropriate  simplifica¬ 
tions.  Despite  the  simplicity  of  the  program,  the  results  agree  well 
with  the  experimental  results.  The  graphs  No.  16  and  No.  17  show 
that  the  constraints  are  well  satisfied  for  both  the  deploy  and  swing 
activities.  Graph  No.  15  depicts  the  gradient  histories  for  the  two 
control  moments  in  the  respective  portions.  In  numerical  computation, 
a  commonly  used  stopping  rule  is  that  the  integral  of  the  square  of  the 
gradients  be  less  than  some  small  positive  number,  i.  e. 


<  e  ;  e  small 


§i  =  du~~  ^°r  swing  and  deploy 

i 

The  condition  that  the  integral  approaches  zero  implies  that  the 
gradient  histories  become  arbitrary  close  to  zero.  From  graph  No. 
15,  the  gradients  are  certainly  not  zero.  These  are  the  best  values 
obtained  from  the  gradient  algorithm  as  the  gradient  trajectories  tend 
to  fluctuate  without  evidence  of  further  improvement  in  the  value  of 
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the  performance  criterion.  The  greatest  deviation  from  zero  occurs 
at  initial  deploy,  when  the  ground  reactions  are  at  their  maximum. 
Judging  from  the  quality  of  our  results,  the  gradient  histories  are 
acceptable  though  not  as  small  as  one  would  have  hoped.  From  our 
computational  experience,  the  cost  improvement  may  be  slight  in 
iterative  calculations  while  the  gradient  histories  may  change  by  a 
wide  margin.  Thus,  the  numerical  limitations  plus  the  large  non¬ 
linearity  of  the  cost  function  account  for  the  gradient  fluctuation  and 
difficulty  in  convergence.  In  addition,  all  our  computations  were 
carried  out  in  single  precision  on  the  IBM  360/65  which  offers 
between  6  and  7  significant  decimal  places.  Double  precision  could 
improve  the  convergence  significantly  if  round  off  errors  are 
accumulating  in  the  computations. 

d.  Equivalence  of  Kinematic  and  Dynamic  Optimality. 

An  important  result  of  the  numerical  simulation  is  the  concept 
of  equivalence  of  kinematic  and  dynamic  optimality.  In  the  absence 
of  the  ground  reactions  and  ankle  moment  (i.  e.  M  =  0,  X  -  0,  Y  ~  0), 
the  model  is  kinematic  in  that  the  moment  contribution  is  entirely 
inertial  and  gravity.  In  the  graphs  No.  13  and  No.  14,  the  f3  =  0  curves 
indicate  the  optimal  moment  profiles  without  any  ground  reaction 
forces.  This  is  the  kinematic  optimal  solution.  The  computation  is 
next  made  for  small  increments  of  the  factor  (3(0^  (3^1)  which 
represent  different  percentages  of  ground  reactions  acting  on  the  foot 
in  stance.  The  graphs  illustrate  the  results  for  the  fractions  (3  -  0. 

0.  55,  0.8  and  1.  0  (the  actual  dynamic  case).  For  each  (3,  the  starting 
control  moments  are  shown  by  the  dotted  lines  obtained  from  a 
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previous  (3*  ((3*  <(3)  curve  by  adding  corrections  due  to  incremental 
ground  reactions.  The  optimal  solutions  are  represented  by  the  solid 
lines.  These  solutions,  which  contain  the  effect  of  ground  reactions, 
are  the  dynamic  optimal  solutions.  On  the  graphs,  the  initial  and 
optimal  curves  are  very  close  together.  These  curves  often  differ  in 
the  first  decimal  place  in  the  computer  output.  Thus,  to  obtain  the 
optimal  moment  profiles  in  the  presence  of  ground  reactions,  one  may 
first  compute  the  optimal  solution  with  no  force  present  and  then  add 
in  the  contribution  from  ground  reactions  and  ankle  moments.  The 
kinematic  solution  ((3  -  0)  is  much  easier  to  obtain  than  the  dynamic 
case  ((3=1)  because  of  the  difficulty  in  handling  numerically  the  large 
forces.  The  following  table  shows  the  (3-series  of  computations. 

Note  that  J  increases  as  the  fraction  of  ground  reactions  is  increased. 
The  higher  values  in  the  "norm"  of  the  gradient  are  probably  due  to 
severe  non-linearity  of  the  cost  surface. 
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e •  Feasible  Initial  and  Terminal  Regions. 

The  starting  and  terminal  values  used  in  the  simulation  are 
chosen  from  experim ental  curves.  In  actual  locomotion,  these 
values  ctr^  seldom  duplicated  either  due  to  parameter  sensitivity  or 
environmental  changes.  Whatever  the  initial  configuration  may  be, 
the  deploy  portion  is  almost  characterized  by  the  same  circular  arc 
equality  constraint.  Thus,  it  would  appear  to  be  very  valuable  to 
obtain  the  initial  and  terminal  feasible  regions  for  adaptation  of  control 
programs.  For  our  purposes,  the  feasible  initial  region  is  defined  as 
the  set  of  values  which  satisfy  the  constraint  and  its  time  derivatives 

at  the  initial  time,  i.e.  Vx  s.t.  S^(x  ;  t  -  0)  =  0  :  S^(x  : t  -  0)  =  0: 

7  o  o7  7  o7  7 

*S^(xoJt=  0)  =  0,  etc.  Results  of  the  calculation  are  shown  in  graphs 
No.  18  to  No.  21.  These  are  obtained  by  varying  the  x^(0)  angle  by 
small  increments  from  the  chosen  value  and  computing  the  corresponding 
x^(0)  through  iteration.  The  striking  fact  is  that  although  the  equation 
is  highly  non-linear,  the  set  of  feasible  values  forms  a  straight  line  as 
shown  in  No.  18.  In  No.  19,  the  velocities  are  also  linearly  related. 

This  is  very  interesting  for  parameter  variation  design  considerations. 
For  a  given  set  of  starting  values,  the  feasible  terminal  values  are 
shown  in  No.  20  and  No.  21. 
f.  Design  Considerations. 

From  the  foregoing  discussions  a.  through  e.  on  the  simulation 
results,  useful  idea^  are  evolved  for  practical  design  and  quantitative 
study  via  the  optimal  programming  approach.  Fig.  19  shows  a 
schematic  realization  of  a  walking  program  by  means  of  multi-arc 
programming.  The  locomotion  behavior  for  a  double -step  is  first 


FIG.  19  REALIZATION  OF  LOCOMOTION  PROGRAM  VIA 


MULTI-ARC  OPTIMAL  PROGRAMMING. 
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decomposed  into  pertinent  sub-modes  or  configurations.  Information 
on  the  feasible  initial  and  terminal  conditions,  temporal  patterns, 
step-length  and  speed  parameters  forms  a  gross  characterization. 

From  this,  one  can  then  generate  the  kinematic  behavior  of  the  gait  - 
i.  e.  displacement,  velocity  and  acceleration  trajectories  for  the 
various  angles.  The  motion  of  the  upper  part  of  the  body  is  likewise 
investigated.  The  trunk  motion  together  with  the  kinematic  information 
specify  the  ground  reactions.  This  leads  to  the  determination  of 
effective  moments,  reactions  at  the  joints  and  other  dynamic  quantities. 
To  facilitate  adaptation  of  the  walking  program,  feedback  of  information 
at  various  levels  is  indicated.  The  linearity  of  the  initial  feasible 
region  as  discussed  in  e.  would  certainly  be  of  value  in  adaptive  design. 
Such  an  optimal  design  procedure  takes  into  account  the  dynamic  and 
kinematic  details  of  the  gait.  This  appears  to  be  much  more  realistic 
than  the  algorithmic  approach  taken  by  McGhee,  Frank  and  Tomovic  in 
their  investigations  [45-48]. 
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ANGULAR  DISPLACEMENT  OF  HIP  AND  KNEE 
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-25. 0L-  ANGULAR  ACCELERATION  OF  HIP  AND  KNEE 


ANGULAR  DISPLACEMENT  (RADIANS) 
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GRAPH  *4 


ROTATION  OF  THE  FOOT  IN  DEPLOY  PORTION 


ANGULAR  ACCELERATIONS  (RAD/SEC2) 
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1.2  — 2.0  SEC'-  SWING 


PHASE  PLANE"  DIAGRAM  OF  HIP  AND  KNEE 
ANGLES  IN  A  DOUBLE-STEP 
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HORIZONTAL  DISPLACEMENT  OF  THE  JOINTS 
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VERTICAL  DISPLACEMENT  OF  THE  JOINTS 
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EFFECTIVE  MOMENTS  AT  THE  JOINTS 
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COMPOSITION  OF  HIP  MOMENT 
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HIP  MOMENT  FOR  VARIOUS  PORTIONS  (/3) 
OF  GROUND  REACTIONS  AT  DEPLOY 
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STATE  EQUALITY  CONSTRAINT  IN  DEPLOY  PORTION 


103 


(j 
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GRAPH  *18 

THE  SET  OF  INITIAL  ANGLES  (x,(0),x2(0)) 
WHICH  SATISFIES  THE  ENTRY  CONDITION 


INITIAL  KNEE  VELOCITY  x2(0)=x4{0)  IN  RAD/SEC 
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GRAPH  *19  . 

THE  SET  OF  INITIAL  VELOCITIES  WHICH  SATISFIES  Sd|  =  0 
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KNEE  ANGLE  x,(ti)  IN  RADIANS 
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HIP  ANGLE  IN  RADIANS  x,(t'f) 


TERMINAL  KNEE  VELOCITY  (tj)  =  X2(t|)  IN  RAD/EEC 
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GRAPH  *21 
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F .  Conclusions 

A  new  quantitative  study  of  the  human  locomotion  problem  has 
been  developed  using  the  methods  of  optimal  programming  and  control 
theory.  The  multi-arc  formulation  is  based  on  the  qualitative 
properties  of  the  biped  gait.  Results  of  the  numerical  experiment 
indicate  that  the  theory  is  highly  relevant  and  agrees  well  with  experi¬ 
mental  results  to-date.  Because  of  the  flexibility  of  the  method,  the 
theory  can  certainly  be  used  to  study  various  types  of  gait,  normal  or 
pathological.  The  quantitative  method  should  find  importance  in 
prosthetic  design.  Eventually,  it  is  hoped,  through  systematic  studies 
like  the  one  described,  locomotion  by  programmed  stimulation  will  be 


realized. 
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In  this  paper,  a  new  approach  to  the  study  of  "programmed"  human  locomotion  is  made 
using  the  methods  of  optimal  programming  and  modern  control  theory. 

Development  of  the  problem  structure  relies  closely  on  the  phasic  and  temporal 
characteristics  of  the  biped  gait.  The  result  is  a  multi -arc  programming  problem  with 
three  stages.  Each  stage  involves  appropriate  dynamic  constraints  which  reflect  the 
particular  nature  of  the  phasic  activity.  Joining  of  the  arcs  is  arranged  in  such  a  way 
as  to  maintain  continuity  of  certain  trajectories  as  well  as  repeatability  of  motion.  A 
novelty  of  the  approach  is  that  the  theory  could  be  used  to  study  walking  behavior  under 
different  environmental  conditions,  such  as  walking  up-stairs  or  over  a  hole. 

A  distinct  feature  of  the  present  approach  which  differs  from  other  studies  is  the 
presence  of  a  minimizing  performance  criterion.  Based  on  external  characteristics  of 
muscles  and  certain  assumptions  regarding  normal  locomotion,  a  simple  quadratic 
type  of  criterion  is  proposed.  This  performance  criterion  is  meaningful  in  that  it  is 
shown  to  be  proportional  to  the  mechanical  work  done  during  normal  walking. 

Invoking  the  necessary  conditions  of  optimal  control  theory,  a  multi-point  boundary 
value  problem  is  obtained.  A  penalty  function  technique  is  employed  for  iterative 
numerical  solution.  Using  the  parameters  available  in  the  literature,  simulation 
results  are  obtained  which  agree  well  with  experimental  studies.  Furthermore, 
refined  features  are  obtained  which  are  not  available  in  previous  studies. 
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